Supplementary Information

# install package librarian if needed
if (!("librarian" %in% rownames(installed.packages()))) {
  install.packages("librarian")
}

# load required packages
librarian::shelf(
  tidyverse,
  usmap,
  fs,
  ggpubr,
  ggdist,
  ggrepel,
  faux,
  lme4,
  lmerTest,
  ggeffects,
  binom,
  tictoc,
  ggthemes,
  sessioninfo,
  knitr,
  kableExtra
)

# Source required functions
myFunctions <- c(
  "FUNStormEventsData_filterData"
)

for (f in myFunctions) {
  source(paste0("../functions/", f, ".R"))
}

# Preperations to show states boundaries
poly_states <- plot_usmap(regions = "states")

# Read in data_details_fips
fileName <- "data_details_fips.RDS"
pathName <- "../data/stormData"
filePath <- dir_ls(path = pathName, regexp = paste0(fileName, "$")) %>% last()
data_details_fips <- readRDS(filePath)

1 Carbon Emission Task: Trials

Table 1 shows all trials in the new variant of the Carbon Emission Task (CET). Each row corresponds to one trial. Carbon Diff. and Bonus Diff. index the unique combinations of relative differences in carbon and bonus consequences in a trial. Carbon Level Rand. and Bonus Level Rand. correspond to the random carbon and bonus level used in the construction of a trial. The average carbon level over all trials is 19.85 lbs CO2, to which we add a random deviation with mean = 0 and SD = 0.9925 lbs CO₂. The average bonus level over all trials is 60 cents, to which we add a random deviation with mean = 0 and SD = 3 cents USD. Carbon Option Pro-Self and Bonus Option Pro-Self represent the actual information available to participants related to the option that will lead to the greater bonus payment for participants. Carbon Option Pro-Climate and Bonus Option Pro-Self represent the same information but for the option leading to the greater benefit for the climate (lower carbon emissions).

# Read in the trials data
CETTrials <- read_csv("../data/CET_trials.csv")

# Select columns to display and round to two digits
CETTrials <- CETTrials %>% 
  select(
    carbon_prcnt, bonus_prcnt,
    carbon_level_rand, bonus_level_rand,
    carbon_self, carbon_env,
    bonus_self, bonus_env
  ) %>% 
  mutate(across(-ends_with("_prcnt"), ~ round(.x, 2)))

# Define column names
CETTrials_colnames <- c(
  "Carbon Diff (%)",
  "Bonus Diff (%)",
  "Carbon Level Rand (lbs.CO2)",
  "Bonus Level Rand (cent USD)",
  "Carbon Option Pro-Self (lbs.CO2)",
  "Carbon Option Pro-Climate (lbs.CO2)",
  "Bonus Option Pro-Self (cent USD)",
  "Bonus Option Pro-Climate (cent USD)"
)

# Dispaly table with custom width for certain columns
kable(CETTrials, col.names = CETTrials_colnames) %>% 
  column_spec(column = 1:2, width = "1cm") %>% 
  column_spec(column = 3:8, width = "2cm")

Table 1: Trials in the new variant of the CET

Carbon Diff (%)	Bonus Diff (%)	Carbon Level Rand (lbs.CO2)	Bonus Level Rand (cent USD)	Carbon Option Pro-Self (lbs.CO2)	Carbon Option Pro-Climate (lbs.CO2)	Bonus Option Pro-Self (cent USD)	Bonus Option Pro-Climate (cent USD)
10	10	20.54	60.33	21.56	19.51	63	57
10	15	19.80	63.84	20.79	18.81	69	59
10	20	19.89	57.64	20.89	18.90	63	52
10	50	20.78	58.24	21.82	19.74	73	44
10	100	20.44	62.12	21.46	19.42	93	31
15	10	20.14	59.66	21.66	18.63	63	57
15	15	19.37	56.33	20.83	17.92	61	52
15	20	18.69	63.63	20.09	17.29	70	57
15	50	19.47	63.24	20.93	18.01	79	47
15	100	18.37	64.15	19.75	16.99	96	32
20	10	19.27	60.23	21.20	17.34	63	57
20	15	18.55	63.58	20.40	16.69	68	59
20	20	19.31	63.65	21.24	17.38	70	57
20	50	19.56	60.05	21.52	17.60	75	45
20	100	20.38	63.13	22.42	18.34	95	32
50	10	19.30	57.58	24.13	14.48	60	55
50	15	20.82	55.45	26.02	15.61	60	51
50	20	20.46	59.72	25.58	15.35	66	54
50	50	20.89	57.35	26.12	15.67	72	43
50	100	21.21	60.94	26.52	15.91	91	30
100	10	19.89	60.80	29.83	9.94	64	58
100	15	20.54	61.23	30.81	10.27	66	57
100	20	19.54	65.30	29.30	9.77	72	59
100	50	20.50	60.10	30.75	10.25	75	45
100	100	19.04	63.07	28.57	9.52	95	32

2 Extreme Weather Data

2.1 Purpose & Rationale

As outlined in the Registered Report, we will assess the number of extreme weather episodes recorded in each participant’s county of residence within the 30 days prior to study completion. Regarding the time window during which we plan to conduct the study, we aim for maximizing the likelihood of capturing suitable variability in the exposure to extreme weather episodes with notable geographic variability. To this end, we analyzed records of extreme weather episodes over the last ten years.

2.2 Filter Data

We filter the storm events data for the specific years, months, and extreme weather event types we are interested in. We filter for all years from 2014 to 2023 (as data are not complete for the year 2024 yet), we highlight the month of July, and we focus on those types of extreme weather events that are predicted to increase in frequency and severity due to climate change : Excessive Heat, Drought, Wildfire, Flash Flood, Coastal Flood, Strong Wind, Hail, and Tornado (IPCC 2023).

# Define variables of interest
myYears <- seq(2014, 2023)
myMonths <- c("July")
myEventTypes <- c(
  "Excessive Heat",
  "Drought",
  "Wildfire",
  "Flash Flood",
  "Coastal Flood",
  "Strong Wind",
  "Hail",
  "Tornado"
)

# Call function
out <- FUNStormEventsData_filterData(
  myData = data_details_fips,
  myYears = myYears,
  myMonths = myMonths,
  myEventTypes = myEventTypes
)

2.3 Analysis

p.hist <- out$dataForHist %>% 
  group_by(year) %>% 
  mutate(
    max_nEpisodes = max(nEpisodes),
    yearlyMean_nEpisodes = mean(nEpisodes)
  ) %>% 
  ungroup() %>% 
  mutate(max_month = ifelse(nEpisodes == max_nEpisodes, TRUE, FALSE)) %>% 
  ggplot(aes(
    x = month_name, y = nEpisodes,
    linewidth = max_month,
    fill = month_name %in% myMonths
  )) +
  geom_hline(
    mapping = aes(yintercept = yearlyMean_nEpisodes),
    linetype = "dashed",
    color = "black"
  ) +
  geom_bar(
    stat = "identity",
    color = "black",
    alpha = .7,
    show.legend = FALSE
  ) +
  scale_linewidth_manual(values = c(0.5, 2)) +
  scale_x_discrete(labels = month.abb) +
  scale_fill_manual(
    values = c("darkgrey", "orange"),
  ) +
  labs(
    title = "Number of Extreme Weather Episodes by Month over the Years 2014 to 2023",
    x = "Month",
    y = "Number of Episodes"
  ) +
  theme_bw() +
  theme(
    text = element_text(size = 15),
    plot.title = element_text(hjust = .5),
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)
  ) +
  facet_wrap(~year, ncol = 5)


jpeg(
  file = "../images/histogramSeasonalDistribution.jpeg",
  width = 14, height = 7.5, units = "in", res = 600
)
print(p.hist)
invisible(dev.off())

p.map_bin <- plot_usmap(
  data = out$dataForUsPlot,
  values = "episodes_bin",
  regions = "counties",
  exclude = c("AK", "HI"),
  color = "black",
  linewidth = 0.1
  ) +
  geom_sf(
    data = poly_states[[1]] %>% 
      filter(!(abbr %in% c("AK", "HI"))),
    color = "black",
    fill = NA,
    linewidth = .3
  ) +
  scale_fill_manual(
    name = "Number of Episodes > 0",
    values = c("white", "orange")
  ) +
  labs(
    title = "Extreme Weather Episodes in July over the Years 2014 to 2023"
  ) +
  theme_bw() +
  theme(
    text = element_text(size = 15),
    legend.position = "bottom",
    plot.title = element_text(hjust = .5),
    panel.grid = element_blank(),
    axis.ticks = element_blank(),
    axis.text = element_blank()
  ) +
  facet_wrap(~year, ncol = 5)

jpeg(
  file = "../images/mapGeographicalDistribution_bin.jpeg",
  width = 14, height = 7.5, units = "in", res = 600
)
print(p.map_bin)
invisible(dev.off())

p.histJulyEventTypes <- data_details_fips %>% 
  mutate(fips = state_county_fips) %>% 
  filter(year %in% myYears) %>% 
  filter(month_name %in% myMonths) %>% 
  filter(event_type %in% myEventTypes) %>% 
  distinct(event_id, fips, .keep_all = TRUE) %>% 
  group_by(year, event_type) %>% 
  summarise(
    nEvents = n(),
    .groups = 'drop'
  ) %>% 
  ggplot(aes(y = event_type, x = nEvents)) +
  geom_bar(
    stat = 'identity',
    color = "black",
    fill = "darkgrey",
    alpha = .7,
  ) +
  geom_text(
    aes(label = scales::comma(nEvents)),
    hjust = -0.1,
    vjust = 0.5,
    size = 2.5
  ) +
  scale_x_continuous(
    labels = scales::comma_format(), 
    breaks = seq(0, 2.5e3, 1e3)
  ) +
  coord_cartesian(xlim = c(0, 2.5e3)) +
  labs(
    title = "Number of Extreme Weather Events in July by Event Type over the Years 2014 to 2023",
    x = "Number of Events",
    y = "Event Type"
  ) +
  theme_bw() +
  theme(
    text = element_text(size = 15),
    plot.title = element_text(hjust = .5)
  ) +
  facet_wrap(~year, ncol = 5)

jpeg(
  file = "../images/histogramJulyEventTypes.jpeg",
  width = 14, height = 7.5, units = "in", res = 600
)
print(p.histJulyEventTypes)
invisible(dev.off())

additionalStats <- data_details_fips %>% 
  filter(year %in% myYears) %>% 
  filter(month_name %in% myMonths) %>% 
  filter(event_type %in% myEventTypes) %>% 
  distinct(episode_id, state_county_fips, .keep_all = TRUE) %>% 
  group_by(year) %>% 
  summarise(
    nEpisodes = n(), 
    .groups = 'drop'
  ) %>% 
  summarise(
    meanNEpisodes = round(mean(nEpisodes), 2),
    sdNEpisodes = round(sd(nEpisodes), 2)
  )

meanNEpisodesJuly <- 
  additionalStats %>% 
  pull(meanNEpisodes) %>% 
  scales::comma()
sdNEpisodesJuly <- 
  additionalStats %>% 
  pull(sdNEpisodes) %>% 
  scales::comma()

Our analysis spanning the last ten years evidenced high numbers of occurrences of extreme weather episodes in every single year (Figure 1) as well as considerable geographical variability (Figure 2). Analyzing the seasonal distribution of extreme weather episodes, Figure 1 shows that July consistently shows a high number of extreme weather episodes over the last ten years, with an average occurrence of 1,777 episodes in each year (SD = 450). Note that numbers represent counts of distinct episode-county combinations to represent how many counties were affected. Table 2 shows the average number of extreme weather events by event type in July over the last ten years, demonstrating that July typically witnesses a multitude of different extreme weather event types. Additionally, Figure 2 indicates that within the month of July, these extreme weather episodes also display a high geographical variability.

Figure 1: Histograms showing the number of extreme weather episodes by month from 2014 to 2023. The dashed horizontal line indicates the mean number of extreme weather episodes in each year. The thick-bordered bar marks the month with the most extreme weather events each year. The orange bar represents July. July had the most extreme weather events in 4 out of 10 years, and in another 4 years, it was right before or after the peak month. Only episodes that included at least one of the following event types were considered: excessive heat, drought, wildfire, flash flood, coastal flood, strong wind, hail, tornado.

data_details_fips %>% 
  mutate(fips = state_county_fips) %>% 
  filter(year %in% myYears) %>% 
  filter(month_name %in% myMonths) %>% 
  filter(event_type %in% myEventTypes) %>% 
  distinct(event_id, fips, .keep_all = TRUE) %>% 
  group_by(year, event_type) %>% 
  summarise(
    nEvents = n(),
    .groups = 'drop'
  ) %>% 
  group_by(event_type) %>% 
  summarise(
    meanNEvents = format(round(mean(nEvents), 2), nsmall = 2),
    sdNEvents = format(round(sd(nEvents, 2)), nsmall = 2),
    .groups = 'drop'
  ) %>% 
  knitr::kable(col.names = c(
    "Event Type", "Mean Number of Events", "SD"
  ), align = c("l", "r", "r"))

Table 2: Mean Number of Extreme Weather Events in July over the Years 2014 to 2023 by Event Type.

Event Type	Mean Number of Events	SD
Coastal Flood	5.43	3.00
Drought	272.70	254.00
Excessive Heat	429.50	347.00
Flash Flood	778.80	241.00
Hail	1249.70	383.00
Strong Wind	8.60	6.00
Tornado	107.60	20.00
Wildfire	106.90	34.00

Figure 2: Maps displaying the geographical distribution of the occurrence of at least one extreme weather episode in July over the years 2014 to 2023. Only episodes that included at least one of the following event types were considered: excessive heat, drought, wildfire, flash flood, coastal flood, strong wind, hail, tornado.

dataForPlot <- out$dataForUsPlot %>% 
  mutate(nEpisodes_withNA = ifelse(nEpisodes == 0, NA_integer_, nEpisodes))

p.map_cont <- plot_usmap(
  data = dataForPlot,
  values = "nEpisodes_withNA",
  regions = "counties",
  exclude = c("AK", "HI"),
  color = "black",
  linewidth = 0.1
  ) +
  geom_sf(
    data = poly_states[[1]] %>% 
      filter(!(abbr %in% c("AK", "HI"))),
    color = "black",
    fill = NA,
    linewidth = .3
  ) +
  scale_fill_binned(
    name = "Number of Episodes",
    n.breaks = 10,
    type = "viridis",
    na.value = "white"
  ) +
  labs(
    title = "Extreme Weather Episodes in July over the Years 2014 to 2023"
  ) +
  theme_bw() +
  theme(
    text = element_text(size = 15),
    legend.position = "bottom",
    plot.title = element_text(hjust = .5),
    panel.grid = element_blank(),
    axis.ticks = element_blank(),
    axis.text = element_blank()
  ) +
  facet_wrap(~year, ncol = 5)

jpeg(
  file = "../images/mapGeographicalDistribution_cont.jpeg",
  width = 14, height = 7.5, units = "in", res = 600
)
print(p.map_cont)
invisible(dev.off())

p.hist_count <- out$dataForUsPlot %>% 
  group_by(year, nEpisodes) %>% 
  summarise(
    count = n(),
    prcnt = count / n_distinct(out$dataForUsPlot$fips)
  ) %>% 
  ggplot(aes(x = nEpisodes, y = prcnt)) +
  geom_bar(stat = "identity", color = "black", fill = "darkgrey") +
  scale_y_continuous(labels = scales::label_percent()) +
  labs(
    x = "Number of Episodes",
    y = "Proportion of Counties"
  ) +
  theme_bw() +
  labs(
    title = "Extreme Weather Episodes in July over the Years 2014 to 2023"
  ) +
  theme(
    text = element_text(size = 15),
    legend.position = "bottom",
    plot.title = element_text(hjust = .5)
  ) +
  facet_wrap(~year, ncol = 5)

jpeg(
  file = "../images/frequencyDistribution_cont.jpeg",
  width = 14, height = 7.5, units = "in", res = 600
)
print(p.hist_count)
invisible(dev.off())

# Calcualte some proportions for display in text
props2023 <- out$dataForUsPlot %>% 
  filter(year == 2023) %>% 
  count(episodes_bin) %>% 
  mutate(
    freq = n/sum(n),
    freq_prcnt = paste0(format(round(freq*100, 2), nsmall = 2), "%")
  )

While Figure 2 visualizes the occurrence of at least one extreme weather episode in July for each county and year (binary variable), Figure 4 displays the actual number of such episodes (continuous). The vast majority of counties were exposed to few episodes, indicating that most of the variability is due to whether an extreme weather episode occurred at all or not. This is further supported by Figure 3 showing histograms for the number of extreme weather episodes in July over the past ten years. Most counties reported either zero or one extreme weather episode in July, and the ratio of counties experiencing no episodes to counties experiencing at least one episode seems to gradually approach 1:1. In July 2023, for instance, this ratio reached 1.02, with 50.43% of counties being exposed to zero and 49.57% of counties being exposed to at least one extreme weather episode.

Figure 3: Maps displaying the geographical distribution of the raw number of extreme weather episodes in July over the years 2014 to 2023. The color palette indicates numbers greater than zero, and white represent a count of zero episodes.

Figure 4: Histograms displaying the distribution of the raw number of extreme weather episodes in July over the years 2014 to 2023. For each number of episodes on the x-axis, the y-axis shows the proportion of counties that recorded this number of episodes.

Finally, as reported in the analysis plan and the design table, we plan to run a set of additional analyses regarding hypotheses H₂ and H₃, in which we will test the sensitivity of results to the time period prior to study completion used to assess extreme weather exposure. Regarding H₂, we will estimate the two-way interaction effect of political affiliation and extreme weather exposure on ΔDuration for different time periods from 30 days to 360 days in increments of 30 days. Similarly for H₃, we will estimate the three-way interaction effect of political affiliation, extreme weather exposure, and attribution of extreme weather events to climate change on ΔDuration for the same time periods. We will visualize results of these additional analyses by plotting the two-way (or three-way) interaction regression coefficients as points surrounded by their 95%-CI on the y-axis and the 12 time periods on the x-axis, as displayed in Figure 5 with simulated data. Based on previous research (Konisky, Hughes, and Kaylor 2016), we expect that the estimated effects will decay as the number of days prior to study completion used to assess the occurrence of extreme weather episodes increases.

set.seed(123)
p.sensitivity <- tibble(
  Days = seq(30, 360, 30),
  Coefficient = accumulate(1:11, ~ .x * .7, .init = 0.1143),
  Error = rnorm(12, .07, 0.005),
  CI_high = Coefficient + .5 * Error,
  CI_low = Coefficient - .5 * Error
) %>% 
  ggplot(aes(x = Days, y = Coefficient, color = CI_low < 0)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_errorbar(aes(ymin = CI_low, ymax = CI_high), width = 3) +
  geom_point(
    shape = "circle filled",
    fill = "white",
    size = 3,
    stroke = 1
  ) +
  scale_color_manual(values = c("black", "grey")) +
  scale_x_continuous(breaks = seq(30, 360, 30)) +
  labs(
    x = "Days used to assess occurence of extreme weather episodes",
    y = "Regression coefficient\n(surrounded by 95%-CI)"
  ) +
  theme_bw() +
  theme(
    panel.grid.minor = element_blank(),
    legend.position = "None"
  )
set.seed(NULL)

jpeg(
  file = "../images/sensitivityAnalyses_simulation.jpeg",
  width = 7, height = 5, units = "in", res = 600
)
print(p.sensitivity)
invisible(dev.off())

Figure 5: Simulated regression coefficients of interaction effects for different number of days prior to study completion used to assess the occurrence of extreme weather episodes. Point estimates are surrounded by their 95% confidence intervals. The dashed line represents the absence of an interaction effect (regression coefficient of zero). Significance of regression coefficients is color-coded, with black indicating regression coefficients significantly different from zero, and grey indicating no significant difference from zero.

2.4 Conclusion

Our analyses indicate that July consistently shows a high number of extreme whether episodes with notable geographic variability (Figure 1 and Figure 2). Therefore, to maximize the likelihood of capturing suitable variability in exposure to extreme weather episodes, we plan to conduct our study at the beginning of August, ensuring that the 30-day period prior to study completion falls within July. Moreover, the main source of variability in exposure to extreme weather episodes in July is due to whether at least one episode occurred or not (Figure 4 and Figure 4). Thus, our main analyses will focus on whether a participant was exposed to at least one extreme weather episode in the 30 days prior to study completion, treated as a binary variable. In additional analyses, we will test the sensitivity of our results to different time periods used to assess extreme weather exposure prior to study completion.

3 Power Analyses

3.1 Purpose & Rationale

The primary goal of the present analyses is to determine the number of participants required to achieve 95% statistical power to detect a Smallest Effect Size Of Interest (SESOI) for a main effect of political affiliation (Democrat vs. Republican) on attentional information search behavior as assessed by ΔDuration. That is, we primarily aim for conducting a power analysis for sample-size determination, also called a priori power analysis (Giner-Sorolla et al. 2024). To this end, we use power simulations with parameters informed by previous studies to assess the statistical power for different sample sizes. This will allow us to decide on a sample size that will provide high statistical power to detect a true and theoretically relevant main effect of political affiliation.

Given this sample size, our secondary goal is to then assess the statistical power to detect different effect sizes for (1) the two-way interaction between political affiliation and extreme weather exposure and (2) the three-way interaction between political affiliation, extreme weather exposure, and the subjective attribution of extreme weather events to climate change. That is, we conduct power analyses for assessing effect-size sensitivity for these two- and three-way interactions (Giner-Sorolla et al. 2024).

The present report is organized as follows:

1. We describe all important variables of the planned study and how we will assess them.

2. We derive a SESOI for the main effect of political affiliation on ΔDuration.

3. We conduct power simulations for sample-size determination analysis based on this political affiliation main effect SESOI.

4. We conduct power simulations for effect-size sensitivity analyses for a two-way interaction effect between political affiliation and extreme weather exposure.

5. We conduct power simulations for effect-size sensitivity analyses for a three-way interaction effect between political affiliation, extreme weather exposure, and subjective attribution of extreme weather events to climate change.

3.2 Important Variables

3.2.1 ΔDuration

Participants will complete 25 trials in a new variant of the Carbon Emission Task (CET, Berger and Wyss 2021) optimized for online process-tracing using MouselabWEB (Willemsen and Johnson 2019). An example trial of the task is displayed and explained in Figure 6.

Figure 6: An example trial of the new variant of the CET. (A) Participants are presented with two options which are associated with different bonus payment and carbon emission consequences. (B) Participants inspect the exact consequences regarding each attribute of each option by hovering their mouse over the respective boxes. Moving the mouse outside of a box occludes the information again. Note that whether the bonus or carbon attributes are displayed in the top row is randomized across participants but held constant within participants. (C) By visiting all boxes, participants can acquire all available information in a trial (all boxes opened simultaneously for demonstration purposes only). Note that whether the option maximizing the bonus payment for participants is presented in the left or right row is randomized within participants.

In all analyses, the criterion (dependent variable) will be ΔDuration (deltaDuration), which is calculated in each trial as:

\[ \Delta Duration = \frac{(t_{Carbon_{A}} + t_{Carbon_{B}}) - (t_{Bonus_{A}} + t_{Bonus_{B}})}{t_{Carbon_{A}} + t_{Carbon_{B}} + t_{Bonus_{A}} + t_{Bonus_{B}}} \]

with \(t\) representing the summed up dwell time on \(Carbon/Bonus\) boxes in Option \(A/B\). ΔDuration varies between -1 and 1 with:

• a value of -1 indicating that the entire dwell time was spent on gathering Bonus information
• a value of 0 indicating that the dwell time was equally split between gathering Bonus and Carbon information
• a value of 1 indication that the entire dwell time was spent on gathering Carbon information

3.2.2 Political Affiliation

We will assess political affiliation using the following question:

Generally speaking, do you think of yourself as a Democrat, Republican or Independent?

Participants will answer on the following 7-point Likert scale:

[1] Strong Republican, [2] Not strong Republican, [3] Independent, close to Republican, [4] Independent, [5] Independent, close to Democrat, [6] Not strong Democrat, [7] Strong Democrat

Additionally, if participants self-identify as [4] Independent, they will be asked:

You said that you think of yourself as an Independent politically. If you had to identify with one party of the two parties, which one would you choose?

Participants will answer on the following scale:

[1] Republican Party, [2] Democratic Party

Based on these questions, we classify participants as either Republican or Democrat as represented by the variable polAff that can take on the following values:

[rep] Republican Party, [dem] Democratic Party

3.2.3 Extreme Weather Exposure

For each participant, we will assess the number of extreme weather episodes that occurred in the participants’ county of residence in the 30 days prior to study completion. We then create a variable ewe which equals to TRUE if at least one extreme weather episode occurred in the specified time interval, and FALSE otherwise.

3.2.4 Subjective Attribution

We will assess the degree to which participants attribute extreme weather events to climate change (see Ogunbode et al. 2019). Participants will rate their agreement to the following three questions:

1. Extreme weather events are caused in part by climate change.

2. Extreme weather events are a sign that the impacts of climate change are happening now.

3. Extreme weather events show us what we can expect from climate change in the future.

Participants will answer on the following 5-point Likert scale:

[1] Strongly disagree, [2] Somewhat disagree, [3] Neither agree nor disagree, [4] Somewhat agree, [5] Strongly agree

Subjective attribution of EWE to climate change will be operationalised as the mean agreement to these three statements. Note that Ogunbode et al. (2019), who used the same questions, response options, and aggregation, report a mean of 3.67 and a SD of 0.85 for this variable.

3.3 SESOI

As argued in the Registered Report, we hypothesize that (H₁) compared to Republicans, Democrats prioritize searching for and attending to carbon over bonus information during decision-making in the CET. That is, Democrats display higher ΔDuration values compared to Republicans.

While H₁ describes the direction of the expected effect, it does not specify the expected magnitude of the effect. This later point is clarified by asking the question what would be the smallest effect size that researchers would still consider theoretically relevant, i.e., what is the smallest effect size of interest (SESOI)? We argue for such a SESOI on theoretical grounds and based on previous MouselabWEB research.

In mouselabWEB studies, it is standard practice to filter out any information acquisition lasting less than 200 msec because such very short (spurious) acquisitions are very unlikely to be consciously processed (Willemsen and Johnson 2019). Therefore, we derive the SESOI based on the consideration that for an effect to be meaningful, the increase in ΔDuration from a Republican to a Democrat should be due to an increase of the time spent gathering carbon (relative to bonus) information of at least 200 msec. We need to translate these considerations into the metric of ΔDuration.

Suppose that Republicans on average spend \(t_{Carbon}\) msec on gathering carbon information and \(t_{Bonus}\) msec on gathering bonus information in each trial, with a total time of gathering any information \(t_{Total} = t_{Carbon} + t_{Bonus}\). Thus, for Republicans, this results in:

\[ \Delta Duration_{Rep} = \frac{t_{Carbon} - t_{Bonus}}{t_{Total}} \]

Suppose that Democrats have the same \(t_{Total}\) as Republicans, but they spend 200 msec more on gathering carbon information and, correspondingly, 200 msec less on gathering bonus information. Thus, for Democrats, this results in:

\[ \Delta Duration_{Dem} = \frac{(t_{Carbon} + 200) - (t_{Bonus} - 200)}{t_{Total}} \\ = \frac{400 + t_{Carbon} - t_{Bonus}}{t_{Total}} \]

Therefore, the difference in \(\Delta Duration\) for Democrats and Republicans is:

\[ \begin{split} \Delta Duration_{Dem} - \Delta Duration_{Rep} = \frac{400 + t_{Carbon} - t_{Bonus}}{t_{Total}} - \frac{t_{Carbon} - t_{Bonus}}{t_{Total}} \\ = \frac {(400 + t_{Carbon} - t_{Bonus}) - (t_{Carbon} - t_{Bonus})}{t_{Total}} \\ = \frac {400}{t_{Total}} \end{split} \]

Thus, we define our SESOI as:

\[SESOI_{polAff} = \frac{400\ msec}{t_{Total}}\]

As this definition shows, the SESOI depends on our expectation regarding \(t_{Total}\). We form these expectations based on previous research. In our lab, we recently conducted another MouselabWEB study whose setup was very similar to the one described in the Registered Report (Studler et al., in preparation). In short, student participants completed a behavioral task in which they searched for and attended to information presented in a 2-by-2 decision matrix adapted from Reeck, Wall, and Johnson (2017). The average time participants spent on acquiring information in each trial (\(t_{Total}\)) was 3.3 seconds. In the planned study we will assess a representative US sample. That is, the average age of our sample will be higher than in a typical student sample. As processing speed is known to decline with age (Salthouse 2000), we assume a slightly higher total information acquisition time in our sample of \(t_{Total} = 3500\ msec\). Therefore, we define our SESOI as:

\[ SESOI_{polAff} = \frac{400\ msec}{3500\ msec} = 0.1143 \]

To account for uncertainty in this estimation, we also provide sample-size determination analysis results based on \(t_{Total} = 4000\ msec\) and \(t_{Total} = 3000\ msec\). This results in a high and low estimate of SESOI:

\[ \begin{split} SESOI_{polAff_{high}} = \frac{400\ msec}{3000\ msec} = 0.1333 \\ SESOI_{polAff_{low}} = \frac{400\ msec}{4000\ msec} = 0.1000 \end{split} \]

3.4 Main Effect

As argued in the Registered Report, we hypothesize:

(H₁) Compared to Republicans, Democrats prioritize searching for and attending to carbon over bonus information during decision-making in the Carbon Emission Task. That is, Democrats display higher ΔDuration values compared to Republicans.

# Create a data frame with predicted true effects

# Smallest Effect Size Of Interest (SESOI)
SESOI <- 0.4/3.5

# Betas
beta_p <- SESOI

# Predicted "true" effects
polAff_trueEffects <- expand_grid(
  polAff = factor(c("rep", "dem"), levels = c("rep", "dem"))
) %>% 
  add_contrast("polAff", contrast = "anova", colnames = "X_p") %>% 
  mutate(
    trueDeltaDuration =
      0 +         # intercept
      X_p * beta_p # main effect polAff
  )

3.4.1 Data Simulation Function

We first define a function that simulates data for the main effect of political affiliation on ΔDuration: FUN_sim. The function will simulate data according to the following model, expressed in lme4-lingo:

deltaDuration ~ polAff + (1|subj) + (1|trial)

The function FUN_sim takes, among others, the following important arguments:

• n_subj: Number of subjects. Chaning this parameter allows us to assess statistical power for different sample sizes.

• beta_0: Fixed intercept (grand mean). We assume that the average participant (irrespective of political affiliation or any other predictor variable) spends about the same time on searching for and attending to carbon as to bonus information in the CET. Therefore, we set beta_0 to zero. The effects of predictor variables will be modeled as deviations from this grand mean.

• beta_p: Fixed effect of political affiliation. This value represents the average difference in ΔDuration between Democrats and Republicans (Democrats - Republicans). As discussed above, we set this value to 0.1143 by default, and we assess the effect of changing this variable on statistical power.

• subj_0: By-subject random intercept SD. We simulate that a participants’ deviations from the grand mean for ΔDuration follows a normal distribution with a mean of 0 and a standard deviation of subj_0 = 0.29. We base our default value on a study by Reeck, Wall, and Johnson (2017). They investigated whether variability in information search behavior is driven predominantly by differences in the features of a choice (i.e., in our case: the relative differences between carbon and bonus outcomes in options A and B) or by individual differences. To this end, they predicted information acquisition using a intercept-only model that included random intercepts for subjects and items. They estimated the random intercept of subjects to be 0.29.

• trial_0: By-trial random intercept SD. We simulate that a items’ deviations from the grand mean for ΔDuration follows a normal distribution with a mean of 0 and a standard deviation of trial_0 = 0.04. Again, we base our default value on the study by Reeck, Wall, and Johnson (2017), who estimated a random intercept for trials of 0.04.

• sigma: Trial-level noise (error) SD. We model the error SD to be of the same size as the sum of the random intercept SDs = 0.29 + 0.04 = 0.33. We also report simulation results for an error SD that is twice the size of the random intercept SDs, i.e., 0.66.

The function FUN_sim is defined below (code visible in html output only):

# define data simulation function
FUN_sim <- function(
  n_subj       =        1000, # number of subjects
  n_subj_prop  =   c(.5, .5), # proportion of republican and democrat subjects
  n_trial      =          25, # number of trials
  beta_0       =           0, # intercept (grand mean) for deltaDuration
  beta_p       =     0.4/3.5, # effect of political affiliation on deltaDuration
  subj_0       =         .29, # by-subject random intercept sd for dt carbon
  trial_0      =         .04, # by-trial random intercept sd
  sigma        = 1*(.29+.04), # residual (error) sd
  
  truncNums    =        TRUE, # should impossible deltaDuration values be truncuated?
  setSeed      =        NULL  # seed number to achieve reproducible results. Set to NULL for simulations!
) {
  
  # set seed to achieve reproducible results for demonstration purposes
  set.seed(setSeed)
  
  # simulate data for dwell time on carbon information
  dataSim <- 
    # add random factor subject
    add_random(subj = n_subj) %>% 
    # add random factor trial
    add_random(trial = n_trial) %>% 
    # add between-subject factor political affiliation (with anova contrast)
    add_between("subj", polAff = c("rep", "dem"), .prob = n_subj_prop*n_subj, .shuffle = FALSE) %>% 
    add_contrast("polAff", colnames = "X_p", contrast = "anova") %>% 
    # add by-subject random intercept
    add_ranef("subj", S_0 = subj_0) %>% 
    # add by-trial random intercept
    add_ranef("trial", T_0 = trial_0) %>% 
    # add error term
    add_ranef(e_st = sigma) %>% 
    # add response values
    mutate(
      # add together fixed and random effects for each effect
      B_0 = beta_0 + S_0 + T_0,
      B_p = beta_p,
      # calculate dv by adding each effect term multiplied by the relevant
      # effect-coded factors and adding the error term
      deltaDuration = B_0 + (B_p * X_p) + e_st
    )
  
  # unset seed
  set.seed(NULL)
  
  # truncuate impossible deltaDurations
  if(truncNums) {
    dataSim <- dataSim %>% 
      mutate(deltaDuration = if_else(deltaDuration < -1, -1,
        if_else(deltaDuration > 1, 1, deltaDuration)))
  }
  
  # run a linear mixed effects model and check summary
  mod <- lmer(
    deltaDuration ~ polAff + (1 | subj) + (1 | trial),
    data = dataSim,
  )
  mod.sum <- summary(mod)
  
  # get results in tidy format
  mod.broom <- broom.mixed::tidy(mod)
  
  return(list(
    dataSim = dataSim,
    modelLmer = mod,
    modelResults = mod.broom
  ))
  
}

We call the function once and extract the results of this single simulation (code visible in html output only):

# Call function
out <- FUN_sim(
  n_subj       =        1000, # number of subjects
  n_subj_prop  =   c(.5, .5), # proportion of republican and democrat subjects
  n_trial      =          25, # number of trials
  beta_0       =           0, # intercept (grand mean) for deltaDuration
  beta_p       =     0.4/3.5, # effect of political affiliation on deltaDuration
  subj_0       =         .29, # by-subject random intercept sd for dt carbon
  trial_0      =         .04, # by-trial random intercept sd
  sigma        = 1*(.29+.04), # residual (error) sd
  
  truncNums    =        TRUE, # should impossible deltaDuration values be truncuated?
  setSeed      =        1234  # seed number to achieve reproducible results. Set to NULL for simulations!
)

# Get results table
resultsTable <- out$modelResults %>% 
  select(-c(std.error, statistic, df)) %>% 
  mutate(across(where(is_double), ~ round(.x, 4))) %>% 
  knitr::kable()
formulaUsedForFit <- paste(as.character(formula(out$modelLmer))[c(2,1,3)], collapse = " ")

# Create plot
p.demo.mainEffect <-  out$dataSim %>% 
  ggplot(aes(x = polAff, y = deltaDuration, color = polAff)) +
  geom_hline(yintercept = 0) +
  geom_violin(alpha = 0.3) +
  geom_point(
    data = polAff_trueEffects,
    mapping = aes(x = polAff, y = trueDeltaDuration),
    shape = "circle open",
    size = 3.5,
    stroke = 2,
    color = "black"
  ) +
  stat_summary(
    fun = mean,
    fun.min = \(x){mean(x) - sd(x)},
    fun.max = \(x){mean(x) + sd(x)},
    position = position_dodge(width = .9)
  ) +
  ggrepel::geom_label_repel(
    data = polAff_trueEffects,
    mapping = aes(x = polAff, y = trueDeltaDuration, label = round(trueDeltaDuration, 4)),
    color = "black",
    box.padding = 1
  ) +
  scale_color_manual(values = c("red", "dodgerblue")) +
  scale_y_continuous(breaks = seq(-1, 1, .2)) +
  labs(title = "Demo Output of One Simulation for Main Effect") +
  theme_bw() +
  theme(legend.position = "none")

jpeg(
  file = "../images/pa_demoMainEffect.jpeg",
  width = 9.1, height = 6.5, units = "in", res = 600
)
print(p.demo.mainEffect)
invisible(dev.off())

Figure 7 visualizes the results of this single simulation and Table 3 summarises the statistical results of fitting the actual model used in data generation to the simulated data.

Figure 7: Visual representation of results of one simulation created using FUN_sim. Violin plots display the full distribution of the data. Points and surrounding lines indicate the mean ± 1 SD. The black horizontal line displays the true sample mean and the black open circles indicate the true means for each cell.

resultsTable

Table 3: Statistical results of one simulation created using FUN_sim. Data was fit using deltaDuration ~ polAff + (1 | subj) + (1 | trial).

effect	group	term	estimate	p.value
fixed	NA	(Intercept)	-0.0188	0.1105
fixed	NA	polAff.dem-rep	0.0896	0.0000
ran_pars	subj	sd__(Intercept)	0.2841	NA
ran_pars	trial	sd__(Intercept)	0.0363	NA
ran_pars	Residual	sd__Observation	0.3223	NA

3.4.2 Power Simulation

We now simulate multiple random samples drawn from the same synthetic population with a known true effect of political affiliation. For each random sample, we fit our statistical model to the data. The statistical power to detect a true effect of political affiliation is calculated as the proportion of significant effects out of the total number of simulations. We aim for a statistical power of 95%.

In the following code, the simulations are calculated. We do not recommend executing this code junk as it takes several hours to run (code only shown in the html version).

FUN_sim_pwr <- function(sim, ...){
  out <- FUN_sim(...)
  modelResults <- out$modelResults %>% 
    mutate(sim = sim) %>% 
    relocate(sim)
  return(modelResults)
}

# How many simulations should be run?
n_sims <- 1000

# What are the breaks for number of subjects we would like to calculate power for?
breaks_subj <- seq(200, 1000, 200)

# What are the breaks for SESOI?
breaks_sesoi <- c(0.4/3, 0.4/3.5, 0.4/4)

# What are the breaks for different errer SDs?
breaks_sigma <- c(1:2)*(.29+.04)

res_mainEffect <- tibble()
for (s in seq_along(breaks_sigma)) {
  
  res_sesoi <- tibble()
  for (sesoi in seq_along(breaks_sesoi)) {
    
    res_nSubj <- tibble()
    for (nSubj in seq_along(breaks_subj)) {
      
      # Give feedback regarding which model is simulated
      cat(paste0(
        "Simulation:\n",
        "  sigma = ", round(breaks_sigma[s], 4), "\n",
        "  sesoi = ", round(breaks_sesoi[sesoi], 4), "\n",
        "  nSubject = ", breaks_subj[nSubj], "\n"
      ))
      
      # Start timer
      cat(paste0("Start date time: ", lubridate::now(), "\n"))
      tic()
      
      # Loop over simulations
      pwr <- map_df(
        1:n_sims, 
        FUN_sim_pwr,
        n_subj = breaks_subj[nSubj],
        beta_p = breaks_sesoi[sesoi],
        sigma = breaks_sigma[s]
      )
      
      # Stop timer and calculate elapsed time
      elapsed_time <- toc(quiet = TRUE)
      elapsed_seconds <- elapsed_time$toc - elapsed_time$tic
      elapsed_minutes <- elapsed_seconds / 60
      cat(paste0("End date time: ", lubridate::now(), "\n"))
      cat("Elapsed time: ", elapsed_minutes, " minutes\n\n")
      
      # Add number of subjects to pwr
      pwr <- pwr %>% 
        mutate(
          nSubjects = breaks_subj[nSubj],
          sesoi = breaks_sesoi[sesoi],
          sigma = breaks_sigma[s]
        )
      
      # Add results to the results table
      res_nSubj <- res_nSubj %>%
        rbind(pwr)
    }
    
    # Add results to the results table
    res_sesoi <- res_sesoi %>% 
      rbind(res_nSubj)
  }
  
  # Add results to the results table
  res_mainEffect <- res_mainEffect %>% 
    rbind(res_sesoi)
  
}

res_mainEffect.summary <- res_mainEffect %>% 
  filter(term == "polAff.dem-rep") %>% 
  group_by(sigma, sesoi, nSubjects) %>% 
  summarise(
    power = mean(p.value < 0.05),
    ci.lower = binom.confint(power*n_sims, n_sims, methods = "exact")$lower,
    ci.upper = binom.confint(power*n_sims, n_sims, methods = "exact")$upper,
    .groups = 'drop'
  ) %>% 
  mutate(
    sigma_fact = factor(format(round(sigma, 4), nsmall = 4)),
    sigma_level = match(sigma_fact, levels(sigma_fact)),
    sesoi_fact = factor(format(round(sesoi, 4), nsmall = 4)),
    sesoi_level = match(sesoi_fact, levels(sesoi_fact))
  )

# Save results in a list object
time <- format(Sys.time(), "%Y%m%d_%H%M")
fileName <- paste0("res_mainEffect", "_", time, ".RDS")
saveRDS(
  list(
    res_mainEffect = res_mainEffect,
    res_mainEffect.summary = res_mainEffect.summary
  ),
  file = file.path("../powerSimulationsOutput", fileName)
)

We retrieve pre-run results:

# Load power simulation data
resList_mainEffect <- readRDS(file.path("../powerSimulationsOutput", "res_mainEffect_20240813_2141.RDS"))
resList_mainEffect.summary <- resList_mainEffect$res_mainEffect.summary

# Extract power values for some specific effect sizes at N = 1000
chosenN <- 600
chosenSigma <- c("0.3300", "0.6600")
chosenSESOI <- c("0.1000", "0.1143")
powerValues <- resList_mainEffect.summary %>% 
  filter(sigma_fact %in% chosenSigma) %>% 
  filter(sesoi_fact %in% chosenSESOI) %>% 
  filter(nSubjects %in% chosenN) %>% 
  mutate(power_str = paste0(round(ci.lower*100, 2), "%")) %>% 
  pull(power_str)

# Extract number of simulations
label_nSimulations <- resList_mainEffect$res_mainEffect$sim %>% n_distinct()

Figure 8 displays the distribution of estimated fixed effects across all simulations. The figure shows that the estimated fixed effects are close to the true ones provided as input in the data simulation function, validating that simulations worked as expected.

Figure 8: Distribution of estimated fixed effects resulting from 1000 simulations for the model deltaDuration ~ polAff + (1 | subj) + (1 | trial). Shaded area represent densities, annotated points indicate medians, and thick and thin lines represent 66% and 95% quantiles.

Figure 9 shows results of our sample-size determination analyses. We find that a sample size of 600 provides a statistical power of 95.4% (lower bound of 95%-CI) even under the most conservative assumptions (SESOI = 0.1000, Error SD = 0.6600).

# Get power values of interest
chosenN <- c(800, 1000)
chosenSigma <- c("0.3300", "0.6600")
chosenSESOI <- c("0.1000")
powerValues <- resList_mainEffect.summary %>% 
  filter(sigma_fact %in% chosenSigma) %>% 
  filter(sesoi_fact %in% chosenSESOI) %>% 
  filter(nSubjects %in% chosenN)

# Get lower CI for power values for conservative sigma assumptions
powerValues_sigmaCons <- powerValues %>% 
  filter(sigma_fact == "0.6600")

# Interpolate the effect sizes at which we achieve 99% power
eff_sigmaCons <- approx(powerValues_sigmaCons$ci.lower, powerValues_sigmaCons$nSubjects, xout = 0.99)$y
# Round results for display in text
eff_sigmaCons_txt <- toString(round(eff_sigmaCons))

Moreover, we inspect the lower bounds of the 95%-CI power estimates in Figure 9. Specifically, we focus on the power simulation results for a conservative SESOI of 0.1000 and an error SD of 0.6600. Then, we use interpolation between sample sizes of 800 and 1000 to estimate the number of participants needed to achieve 99% statistical power (lower bounds of 95%-CI power estimates). Results reveal that with a sample size of N = 895, we would achieve 99% statistical power (lower bound of 95%-CI power estimate) to detect a true effect of at least 0.1000.

We optimize the study design to detect a true SESOI for political affiliation. However, we are also interested in two-way and three-way interaction effects, which are known to require greater sample sizes to achieve the same statistical power as for main effects. Moreover, greater sample sizes are more likely to accurately represent target populations with respect to variables like exposure to extreme weather events and subjective attribution of extreme weather events to climate change. Therefore, we opt for a sample size of N = 1000 for further effect-size sensitivity analyses regarding interaction effects.

Figure 9: Power curves for the main effect of polAff. Points represent statistical power surrounded by a 95%-CI based on 1000 simulations with α = 0.05.

3.5 Two-Way Interaction Effect

As argued in the Registered Report, we hypothesize:

(H₂): Exposure to extreme weather events amplifies the polarization in attentional information search processes, increasing the difference in ΔDuration between Republicans and Democrats among individuals with (vs. without) such exposure. In other words, there is a positive two-way interaction effect of extreme weather exposure and political affiliation on ΔDuration.

# Create a data frame with predicted true effects

# Smallest Effect Size Of Interest (SESOI)
SESOI <- 0.4/3.5

# Betas
beta_p <- SESOI
beta_e <- 0
beta_p_e_inx <- SESOI

# Predicted "true" effects
polAff_ewe_trueEffects <- expand_grid(
  polAff = factor(c("rep", "dem"), levels = c("rep", "dem")),
  ewe = factor(c(FALSE, TRUE), levels = c(FALSE, TRUE))
) %>% 
  add_contrast("polAff", contrast = "anova", colnames = "X_p") %>% 
  add_contrast("ewe", contrast = "anova", colnames = "X_e") %>% 
  mutate(
    trueDeltaDuration =
      0 +                        # intercept
      X_p * beta_p +             # main effect polAff
      X_e * beta_e +             # main effect ewe
      X_p * X_e * beta_p_e_inx   # interaction effect polAff*ewe
  )

3.5.1 SESOI for Two-Way Interaction

For deriving a SESOI for the two-way interaction of interest, similar considerations apply as in the case of the SESOI for the main effect of interest. In this section, we make these considerations explicit. We start by noticing that the complete fixed two-way interaction polAff × ewe is modeled as:

\[ \Delta Duration = \beta_{0} + \beta_{1} \cdot polAff + \beta_{2} \cdot ewe + \beta_{3} \cdot (polAff \times ewe) \]

By rearranging terms, one can show that the effect of polAff is given by:

\[ Effect_{polAff} = \beta_{1} + \beta_{3} \cdot ewe \]

Now, let’s calculate this effect for two individuals who differ in their levels of ewe. As ewe is a binary variable, we have two types of individuals: individuals with (ewe = 1) and without (ewe = 0) extreme weather exposure. For an individual without extreme weather exposure, the effect of polAff will be:

\[ Effect_{noEWE} = \beta_{1} + \beta_{3} \cdot 0 = \beta_{1} \]

For an individual with extreme weather exposure, the effect of polAff will be:

\[ Effect_{EWE} = \beta_{1} + \beta_{3} \cdot 1 = \beta_{1} + \beta_{3} \]

The difference in the effect of polAff between these two individuals is given by:

\[ \begin{split} Effect_{EWE} - Effect_{noEWE} = (\beta_{1} + \beta_{3}) - \beta_{1} = \beta_{3} \end{split} \]

We consider this difference as theoretically relevant if it is at least of the same size as the SESOI for the effect of polAff:

\[ SESOI_{polAff \times ewe} = Effect_{EWE} - Effect_{noEWE} = SESOI_{polAff} = 0.1143 \]

3.5.2 Data Simulation Function

We next define a function that simulates data for the two-way interaction effect of political affiliation with extreme weather exposure on ΔDuration: FUN_sim_2wayInt. The function will simulate data according to the following model, expressed in lme4-lingo:

deltaDuration ~ polAff * ewe + (1|subj) + (1|trial)

The function FUN_sim_2wayInt takes, among others, the following important arguments (in addition to the arguments discussed for FUN_sim):

• beta_p: Fixed main effect of political affiliation. Compatible with H₁, we keep this value at the SESOI of 0.1143. That is, we model that the average effect of political affiliation across all participants, irrespective of their extreme weather exposure, is 0.1143.

• beta_e: Fixed main effect of extreme weather exposure. As we are interested in the moderating role of extreme weather exposure, we set this main effect to zero. That is, we assume that the effect of extreme weather exposure is highly dependent on participants’ political affiliation.

• beta_p_e_inx: Fixed two-way interaction effect of political affiliation and extreme weather exposure. We set this initial value to the SESOI derived above, but we investigate how changing this effect size impacts statistical power, as we are conducting effect-size sensitivity analyses for interaction effects.

The function FUN_sim_2wayInt is defined below:

# define data simulation function
FUN_sim_2wayInt <- function(
  n_subj         =        1000, # number of subjects
  n_subj_prop_p  =   c(.5, .5), # proportion of republican and democrat subjects
  n_subj_prop_e  =   c(.5, .5), # proportion of subjects without and with extreme weather exposure
  n_trial        =          25, # number of trials
  beta_0         =           0, # intercept (grand mean) for deltaDuration
  beta_p         =     0.4/3.5, # main effect of political affiliation (polAff)
  beta_e         =           0, # main effect of extreme weather exposure (ewe)
  beta_p_e_inx   =     0.4/3.5, # two-way interaction effect of polAff and ewe
  subj_0         =         .29, # by-subject random intercept sd for dt carbon
  trial_0        =         .04, # by-trial random intercept sd
  sigma          = 1*(.29+.04), # residual (error) sd
  
  truncNums      =        TRUE, # should impossible numbers be truncuated?
  setSeed        =        NULL  # seed number to achieve reproducible results. Set to NULL for simulations!
) {
  
  # set seed to achieve reproducible results for demonstration purposes
  set.seed(setSeed)
  
  # simulate data for dwell time on carbon information
  dataSim <- 
    # add random factor subject
    add_random(subj = n_subj) %>% 
    # add random factor trial
    add_random(trial = n_trial) %>% 
    # add between-subject factor political affiliation (with anova contrast)
    add_between("subj", polAff = c("rep", "dem"), .prob = n_subj_prop_p*n_subj, .shuffle = TRUE) %>% 
    add_contrast("polAff", colnames = "X_p", contrast = "anova") %>% 
    # add between-subject factor extreme weather exposure (with anova contrast)
    add_between("subj", ewe = c(FALSE, TRUE), .prob = n_subj_prop_e*n_subj, .shuffle = TRUE) %>% 
    add_contrast("ewe", colnames = "X_e", contrast = "anova") %>% 
    # add by-subject random intercept
    add_ranef("subj", S_0 = subj_0) %>% 
    # add by-trial random intercept
    add_ranef("trial", T_0 = trial_0) %>% 
    # add error term
    add_ranef(e_st = sigma) %>% 
    # add response values
    mutate(
      # add together fixed and random effects for each effect
      B_0 = beta_0 + S_0 + T_0,
      B_p = beta_p,
      B_e = beta_e,
      B_p_e_inx = beta_p_e_inx,
      # calculate dv by adding each effect term multiplied by the relevant
      # effect-coded factors and adding the error term
      deltaDuration = 
        B_0 + e_st +
        (X_p * B_p) +
        (X_e * B_e) +
        (X_p * X_e * B_p_e_inx)
    )
  
  # truncuate impossible deltaDurations
  if(truncNums) {
    dataSim <- dataSim %>% 
      mutate(deltaDuration = if_else(deltaDuration < -1, -1,
        if_else(deltaDuration > 1, 1, deltaDuration)))
  }
  
  # run a linear mixed effects model and check summary
  mod <- lmer(
    deltaDuration ~ polAff*ewe + (1 | subj) + (1 | trial),
    data = dataSim
  )
  mod.sum <- summary(mod)
  
  # get results in tidy format
  mod.broom <- broom.mixed::tidy(mod)
  
  return(list(
    dataSim = dataSim,
    modelLmer = mod,
    modelResults = mod.broom
  ))
  
}

We call the function once and extract the results of this single simulation:

out <- FUN_sim_2wayInt(
  n_subj         =        1000, # number of subjects
  n_subj_prop_p  =   c(.5, .5), # proportion of republican and democrat subjects
  n_subj_prop_e  =   c(.5, .5), # proportion of subjects without and with extreme weather exposure
  n_trial        =          25, # number of trials
  beta_0         =           0, # intercept (grand mean) for deltaDuration
  beta_p         =     0.4/3.5, # main effect of political affiliation (polAff)
  beta_e         =           0, # main effect of extreme weather exposure (ewe)
  beta_p_e_inx   =     0.4/3.5, # two-way interaction effect of polAff and ewe
  subj_0         =         .29, # by-subject random intercept sd for dt carbon
  trial_0        =         .04, # by-trial random intercept sd
  sigma          = 1*(.29+.04), # residual (error) sd
  
  truncNums      =        TRUE, # should impossible numbers be truncuated?
  setSeed        =        1234  # seed number to achieve reproducible results. Set to NULL for simulations!
)

# Get results table
resultsTable <- out$modelResults %>% 
  select(-c(std.error, statistic, df)) %>% 
  mutate(across(where(is_double), ~ round(.x, 4))) %>% 
  knitr::kable()
formulaUsedForFit <- paste(as.character(formula(out$modelLmer))[c(2,1,3)], collapse = " ")

# Create plot
p.demo.2wayInt <-  out$dataSim %>% 
  ggplot(aes(x = ewe, y = deltaDuration, color = polAff)) +
  geom_hline(yintercept = 0) +
  geom_violin(alpha = 0.3) +
  geom_point(
    data = polAff_ewe_trueEffects,
    mapping = aes(x = ewe, y = trueDeltaDuration, fill = polAff),
    shape = "circle open",
    size = 3.5,
    stroke = 2,
    color = "black",
    position = position_dodge(width = .9),
    show.legend = FALSE
  ) +
  stat_summary(
    fun = mean,
    fun.min = \(x){mean(x) - sd(x)},
    fun.max = \(x){mean(x) + sd(x)},
    position = position_dodge(width = .9)
  ) +
  ggrepel::geom_label_repel(
    data = polAff_ewe_trueEffects,
    mapping = aes(x = ewe, y = trueDeltaDuration, fill = polAff, label = round(trueDeltaDuration, 4)),
    color = "black",
    box.padding = 1,
    position = position_dodge(width = .9),
    show.legend = FALSE
  ) +
  scale_color_manual(values = c("red", "dodgerblue")) +
  scale_fill_manual(values = c("white", "white")) +
  scale_y_continuous(breaks = seq(-1, 1, .2)) +
  labs(title = "Demo Output of One Simulation for Two-Way Interaction") +
  theme_bw()

jpeg(
  file = "../images/pa_demo2wayInt.jpeg",
  width = 9.1, height = 6.3, units = "in", res = 600
)
print(p.demo.2wayInt)
invisible(dev.off())

Figure 10 visualizes the results of this single simulation and Table 4 summarizes the statistical results of fitting the actual model used in data generation to the simulated data.

Figure 10: Visual representation of results of one simulation created using FUN_sim_2wayInt. Violin plots display the full distribution of the data. Points and surrounding lines indicate the mean ± 1 SD. The black horizontal line displays the true sample mean and the black open circles indicate the true means for each cell.

Table 4: Statistical results of one simulation created using FUN_sim_2wayInt. Data was fit using deltaDuration ~ polAff * ewe + (1 | subj) + (1 | trial).

effect	group	term	estimate	p.value
fixed	NA	(Intercept)	-0.0006	0.9625
fixed	NA	polAff.dem-rep	0.1290	0.0000
fixed	NA	ewe.TRUE-FALSE	-0.0139	0.4516
fixed	NA	polAff.dem-rep:ewe.TRUE-FALSE	0.1398	0.0002
ran_pars	subj	sd__(Intercept)	0.2855	NA
ran_pars	trial	sd__(Intercept)	0.0382	NA
ran_pars	Residual	sd__Observation	0.3232	NA

3.5.3 Power Simulation

In the following code, the simulations are calculated. We do not recommend executing this code junk as it takes several hours to run.

FUN_sim_2wayInt_pwr <- function(sim, ...){
  out <- FUN_sim_2wayInt(...)
  modelResults <- out$modelResults %>% 
    mutate(sim = sim) %>% 
    relocate(sim)
  return(modelResults)
}

# How many simulations should be run?
n_sims <- 1000

# What are the breaks for number of subjects we would like to calculate power for?
breaks_subj <- c(900, 950, 1000)

# What are the breaks for SESOI?
breaks_sesoi <- (0.4/3.5)*seq(1, 2, .25)

# What are the breaks for different error SDs?
breaks_sigma <- c((.29+.04), 2*(.29+.04))

res_2wayInt <- tibble()
for (s in seq_along(breaks_sigma)) {
  
  res_sesoi <- tibble()
  for (sesoi in seq_along(breaks_sesoi)) {
    
    res_nSubj <- tibble()
    for (nSubj in seq_along(breaks_subj)) {
      
      # Give feedback regarding which model is simulated
      cat(paste0(
        "Simulation:\n",
        "  sigma = ", round(breaks_sigma[s], 4), "\n",
        "  sesoi = ", round(breaks_sesoi[sesoi], 4), "\n",
        "  nSubject = ", breaks_subj[nSubj], "\n"
      ))
      
      # Start timer
      cat(paste0("Start date time: ", lubridate::now(), "\n"))
      tic()
      
      # Loop over simulations
      pwr <- map_df(
        1:n_sims, 
        FUN_sim_2wayInt_pwr,
        n_subj = breaks_subj[nSubj],
        beta_p_e_inx = breaks_sesoi[sesoi],
        sigma = breaks_sigma[s]
      )
      
      # Stop timer and calculate elapsed time
      elapsed_time <- toc(quiet = TRUE)
      elapsed_seconds <- elapsed_time$toc - elapsed_time$tic
      elapsed_minutes <- elapsed_seconds / 60
      cat(paste0("End date time: ", lubridate::now(), "\n"))
      cat("Elapsed time: ", elapsed_minutes, " minutes\n\n")
      
      # Add number of subjects to pwr
      pwr <- pwr %>% 
        mutate(
          nSubjects = breaks_subj[nSubj],
          sesoi = breaks_sesoi[sesoi],
          sigma = breaks_sigma[s]
        )
      
      # Add results to the results table
      res_nSubj <- res_nSubj %>%
        rbind(pwr)
    }
    
    # Add results to the results table
    res_sesoi <- res_sesoi %>% 
      rbind(res_nSubj)
  }
  
  # Add results to the results table
  res_2wayInt <- res_2wayInt %>% 
    rbind(res_sesoi)
  
}

res_2wayInt.summary <- res_2wayInt %>% 
  filter(term == "polAff.dem-rep:ewe.TRUE-FALSE") %>% 
  group_by(sigma, sesoi, nSubjects) %>% 
  summarise(
    power = mean(p.value < 0.05),
    ci.lower = binom.confint(power*n_sims, n_sims, methods = "exact")$lower,
    ci.upper = binom.confint(power*n_sims, n_sims, methods = "exact")$upper,
    .groups = 'drop'
  ) %>% 
  mutate(
    sigma_fact = factor(format(round(sigma, 4), nsmall = 4)),
    sigma_level = match(sigma_fact, levels(sigma_fact)),
    sesoi_fact = factor(format(round(sesoi, 4), nsmall = 4)),
    sesoi_level = match(sesoi_fact, levels(sesoi_fact))
  )

# Save results in a list object
time <- format(Sys.time(), "%Y%m%d_%H%M")
fileName <- paste0("res_2wayInt", "_", time, ".RDS")
saveRDS(
  list(
    res_2wayInt = res_2wayInt,
    res_2wayInt.summary = res_2wayInt.summary
  ),
  file = file.path("../powerSimulationsOutput", fileName)
)

We retrieve pre-run results:

# Load power simulation data
resList_2wayInt <- readRDS(file.path("../powerSimulationsOutput", "res_2wayInt_20240814_1052.RDS"))
resList_2wayInt.summary <- resList_2wayInt$res_2wayInt.summary

# Extract power values for some specific assumptions
chosenN <- 950
chosenSigma <- c("0.3300", "0.6600")
chosenSESOI <- c("0.1429", "0.1714")
powerValues <- resList_2wayInt.summary %>% 
  filter(sigma_fact %in% chosenSigma) %>% 
  filter(sesoi_fact %in% chosenSESOI) %>% 
  filter(nSubjects %in% chosenN) %>% 
  mutate(power_str = paste0(round(power*100, 2), "%")) %>% 
  pull(power_str)

# Extract number of simulations
label_nSimulations <- resList_2wayInt$res_2wayInt$sim %>% n_distinct()

# Repeat breaks_sesoi
breaks_sesoi <- (0.4/3.5)*seq(1, 2, .25)

Figure 11 displays the distribution of estimated fixed effects across all simulations. The figure shows that the estimated fixed effects are close to the true ones provided as input in the data simulation function, validating that simulations worked as expected.

# Define some filters 
filter_nSubjects <- 1000
filter_sesoi <- unique(resList_2wayInt$res_2wayInt$sesoi)[1]
filter_sigma <- unique(resList_2wayInt$res_2wayInt$sigma)[1]

# Prepare data for plot
fixedEstimates <- resList_2wayInt$res_2wayInt %>% 
  mutate(
    group = ifelse(is.na(group), "", group),
    group_term = str_remove(str_c(group, term, sep = "_"), "^_")
  ) %>% 
  filter(effect == "fixed") %>% 
  filter(nSubjects == filter_nSubjects) %>% 
  filter(sesoi == filter_sesoi) %>% 
  filter(sigma == filter_sigma)
fixedEstimates_medians <- fixedEstimates %>% 
  group_by(group_term) %>% 
  summarise(
    median = median(estimate, na.rm = TRUE),
    median_rounded = format(round(median, 4), nsmall = 4, scientific = FALSE),
    .groups = 'drop'
  )

# Create plot
p.checkSims.2wayInt <-fixedEstimates %>% 
  ggplot(aes(x = estimate, y = group_term)) +
  ggdist::stat_halfeye() +
  ggrepel::geom_label_repel(
    data = fixedEstimates_medians,
    mapping = aes(x = median, y = group_term, label = median_rounded),
    box.padding = .5
  ) +
  labs(
    title = "Distribution of Estimated Fixed Effects",
    subtitle = paste0(
      "SESOI = ", round(filter_sesoi, 4), ", ",
      "σ = ", round(filter_sigma, 4), ", ",
      "N = ", filter_nSubjects, ", ",
      "Number of Simulations = ", label_nSimulations
    ),
    x = "Estimate",
    y = "Term"
  ) +
  theme_bw()

jpeg(
  file = "../images/pa_checkSims2wayInt.jpeg",
  width = 9.1, height = 6.5, units = "in", res = 600
)
print(p.checkSims.2wayInt)
invisible(dev.off())

Figure 11: Distribution of estimated fixed effects resulting from 1000 simulations for the model deltaDuration ~ polAff * ewe + (1 | subj) + (1 | trial). Shaded area represent densities, annotated points indicate medians, and thick and thin lines represent 66% and 95% quantiles.

Figure 12 shows results of our effect-size sensitivity analyses. We plot statistical power (y-axis) for different effect sizes (x-axis), taking into account different assumptions for the error SD (color) and sample size (panel). Regarding the latter, we report results not only for the full sample size we aim for (N = 1000), but also for sample sizes taking into account different participant exclusion-rates due to exclusion criteria defined in the Registered Report.

Figure 12: Power curves for the two-way interaction polAff × ewe. Points represent simulated power surrounded by a 95%-CI based on 1000 simulations with α = 0.05. Note that, in contrast to Figure 9, the x-axsis represents different effect sizes, starting from the defined SESOI, while the panels represent different sample sizes, taking into account participant exclusion-rates of 10% (N = 900), 5% (N = 950), and 0% (N = 1000). Note that estimates are displayed with a slight shift along the x-axis to reduce overlap. Since we only consider effect sizes equal to or greater than the SESOI to be of relevance, power curves were only examined for positive effect sizes, which is why only positive values are shown.

# Get power values of interest
chosenN <- 950
chosenSigma <- c("0.3300", "0.6600")
chosenSESOI <- c("0.1429", "0.1714")
powerValues <- resList_2wayInt.summary %>% 
  filter(sigma_fact %in% chosenSigma) %>% 
  filter(sesoi_fact %in% chosenSESOI) %>% 
  filter(nSubjects %in% chosenN)

# Get lower CI for power values for more liberal and more conservative sigma assumptions
powerValues_sigmaLib <- powerValues %>% 
  filter(sigma_fact == "0.3300")
powerValues_sigmaCons <- powerValues %>% 
  filter(sigma_fact == "0.6600")

# Interpolate the effect sizes at which we achieve 95% power
eff_sigmaLib <- approx(powerValues_sigmaLib$ci.lower, powerValues_sigmaLib$sesoi, xout = 0.95)$y
eff_sigmaCons <- approx(powerValues_sigmaCons$ci.lower, powerValues_sigmaCons$sesoi, xout = 0.95)$y
# Round results for display in text
eff_sigmaLib_txt <- format(round(eff_sigmaLib, 4), nsmall = 4)
eff_sigmaCons_txt <- format(round(eff_sigmaCons, 4), nsmall = 4)

To assess the smallest effect size that can be detected with 95% statistical power, we inspect the lower bounds of the 95%-CI power estimates in Figure 12. Specifically, we focus on the power simulation results for N = 950, which takes into account a participant exclusion-rate of 5%. There, we interpolate between the two point estimates that lie just below and above the 95% power line, i.e., between the power estimates for effect sizes 0.1429 and 0.1714. Assuming an error SD of 0.3300, we achieve 95% statistical power to detect a two-way interaction effect of at least 0.1530. For a more conservative error SD of 0.6600, this smallest detectable effect size is only marginally higher (0.1619).

3.5.3.1 Increased Sample Size

We additionally conduct an effect-size sensitivity analysis for increased sample sizes.

FUN_sim_2wayInt_pwr <- function(sim, ...){
  out <- FUN_sim_2wayInt(...)
  modelResults <- out$modelResults %>% 
    mutate(sim = sim) %>% 
    relocate(sim)
  return(modelResults)
}

# How many simulations should be run?
n_sims <- 1000

# What are the breaks for number of subjects we would like to calculate power for?
breaks_subj <- c(1500, 1750, 2000)

# What are the breaks for SESOI?
breaks_sesoi <- (0.4/3.5)*seq(1, 2, .25)

# What are the breaks for different error SDs?
breaks_sigma <- c((.29+.04), 2*(.29+.04))

res_2wayInt <- tibble()
for (s in seq_along(breaks_sigma)) {
  
  res_sesoi <- tibble()
  for (sesoi in seq_along(breaks_sesoi)) {
    
    res_nSubj <- tibble()
    for (nSubj in seq_along(breaks_subj)) {
      
      # Give feedback regarding which model is simulated
      cat(paste0(
        "Simulation:\n",
        "  sigma = ", round(breaks_sigma[s], 4), "\n",
        "  sesoi = ", round(breaks_sesoi[sesoi], 4), "\n",
        "  nSubject = ", breaks_subj[nSubj], "\n"
      ))
      
      # Start timer
      cat(paste0("Start date time: ", lubridate::now(), "\n"))
      tic()
      
      # Loop over simulations
      pwr <- map_df(
        1:n_sims, 
        FUN_sim_2wayInt_pwr,
        n_subj = breaks_subj[nSubj],
        beta_p_e_inx = breaks_sesoi[sesoi],
        sigma = breaks_sigma[s]
      )
      
      # Stop timer and calculate elapsed time
      elapsed_time <- toc(quiet = TRUE)
      elapsed_seconds <- elapsed_time$toc - elapsed_time$tic
      elapsed_minutes <- elapsed_seconds / 60
      cat(paste0("End date time: ", lubridate::now(), "\n"))
      cat("Elapsed time: ", elapsed_minutes, " minutes\n\n")
      
      # Add number of subjects to pwr
      pwr <- pwr %>% 
        mutate(
          nSubjects = breaks_subj[nSubj],
          sesoi = breaks_sesoi[sesoi],
          sigma = breaks_sigma[s]
        )
      
      # Add results to the results table
      res_nSubj <- res_nSubj %>%
        rbind(pwr)
    }
    
    # Add results to the results table
    res_sesoi <- res_sesoi %>% 
      rbind(res_nSubj)
  }
  
  # Add results to the results table
  res_2wayInt <- res_2wayInt %>% 
    rbind(res_sesoi)
  
}

res_2wayInt.summary <- res_2wayInt %>% 
  filter(term == "polAffX_p:eweX_e") %>% 
  group_by(sigma, sesoi, nSubjects) %>% 
  summarise(
    power = mean(p.value < 0.05),
    ci.lower = binom.confint(power*n_sims, n_sims, methods = "exact")$lower,
    ci.upper = binom.confint(power*n_sims, n_sims, methods = "exact")$upper,
    .groups = 'drop'
  ) %>% 
  mutate(
    sigma_fact = factor(format(round(sigma, 4), nsmall = 4)),
    sigma_level = match(sigma_fact, levels(sigma_fact)),
    sesoi_fact = factor(format(round(sesoi, 4), nsmall = 4)),
    sesoi_level = match(sesoi_fact, levels(sesoi_fact))
  )

# Save results in a list object
time <- format(Sys.time(), "%Y%m%d_%H%M")
fileName <- paste0("res_2wayInt_increasedN", "_", time, ".RDS")
saveRDS(
  list(
    res_2wayInt = res_2wayInt,
    res_2wayInt.summary = res_2wayInt.summary
  ),
  file = file.path("../powerSimulationsOutput", fileName)
)

# Load power simulation data
resList_2wayInt <- readRDS(file.path("../powerSimulationsOutput", "res_2wayInt_increasedN_20250308_0703.RDS"))
resList_2wayInt.summary <- resList_2wayInt$res_2wayInt.summary

# Extract power values for some specific assumptions
chosenN <- 2000
chosenSigma <- c("0.3300", "0.6600")
chosenSESOI <- c("0.1429", "0.1714")
powerValues <- resList_2wayInt.summary %>% 
  filter(sigma_fact %in% chosenSigma) %>% 
  filter(sesoi_fact %in% chosenSESOI) %>% 
  filter(nSubjects %in% chosenN) %>% 
  mutate(power_str = paste0(round(power*100, 2), "%")) %>% 
  pull(power_str)

# Extract number of simulations
label_nSimulations <- resList_2wayInt$res_2wayInt$sim %>% n_distinct()

# Repeat breaks_sesoi
breaks_sesoi <- (0.4/3.5)*seq(1, 2, .25)

Figure 12 shows results of our effect-size sensitivity analyses. We plot statistical power (y-axis) for different effect sizes (x-axis), taking into account different assumptions for the error SD (color) and sample size (panel). Regarding the latter, we report results not only for the full sample size we aim for (N = 1000), but also for sample sizes taking into account different participant exclusion-rates due to exclusion criteria defined in the Registered Report.

Figure 13: Power curves for the two-way interaction polAff × ewe. Points represent simulated power surrounded by a 95%-CI based on 1000 simulations with α = 0.05. Note that, in contrast to Figure 9, the x-axsis represents different effect sizes, starting from the defined SESOI, while the panels represent different sample sizes. Note that estimates are displayed with a slight shift along the x-axis to reduce overlap. Since we only consider effect sizes equal to or greater than the SESOI to be of relevance, power curves were only examined for positive effect sizes, which is why only positive values are shown.

3.6 Three-Way Interaction Effect

As argued in the Registered Report, we hypothesize:

(H₃): The two-way interaction effect of extreme weather exposure and political affiliation is greater for individuals who more strongly attribute extreme weather events to climate change. In other words, there is a positive three-way interaction effect of political affiliation, exposure to extreme weather events, and their attribution to climate change on ΔDuration.

# Smallest Effect Size Of Interest (SESOI)
SESOI <- 0.4/3.5

# Betas
beta_p <- SESOI
beta_e <- 0
beta_s <- 0
beta_p_e_inx <- SESOI
beta_p_s_inx <- 0
beta_e_s_inx <- 0
beta_p_e_s_inx <- SESOI

# Predicted "true" effects
polAff_ewe_subjAttr_trueEffects <- expand_grid(
  polAff = factor(c("rep", "dem"), levels = c("rep", "dem")),
  ewe = factor(c(FALSE, TRUE), levels = c(FALSE, TRUE)),
  subjAttr = factor(c("-SD", "M", "+SD"), levels = c("-SD", "M", "+SD"))
) %>% 
  add_contrast("polAff", contrast = "anova", colnames = "X_p") %>% 
  add_contrast("ewe", contrast = "anova", colnames = "X_e") %>% 
  add_contrast("subjAttr", contrast = "anova", colnames = c("X_s_1", "X_s_2")) %>% 
  mutate(
    trueDeltaDuration =
      0 +                                      # intercept
      X_p * beta_p +                           # main effect polAff
      X_e * beta_e +                           # main effect ewe
      X_s_1 * beta_s +                         # main effect subjAttr, dummy variable for M vs. -SD
      X_s_2 * (2 * beta_s) +                   # main effect subjAttr, dummy variable for +SD vs. -SD
      X_p * X_e * beta_p_e_inx +               # 2-way interaction polAff * ewe
      X_p * X_s_1 * beta_p_s_inx +             # 2-way interaction polAff * subjAttr, dummy variable for M vs. -SD
      X_p * X_s_2 * (2 * beta_p_s_inx) +       # 2-way interaction polAff * subjAttr, dummy variable for +SD vs. -SD
      X_e * X_s_1 * beta_e_s_inx +             # 2-way interaction ewe * subjAttr, dummy variable for M vs. -SD
      X_e * X_s_2 * (2 * beta_e_s_inx) +       # 2-way interaction ewe * subjAttr, dummy variable for +SD vs. -SD
      X_p * X_e * X_s_1 * beta_p_e_s_inx +     # 3-way interaction polAff*ewe*subjAttr, dummy variable for M vs -SD
      X_p * X_e * X_s_2 * (2 * beta_p_e_s_inx) # 3-way interaction polAff*ewe*subjAttr, dummy variable for +SD vs. -SD
  )

3.6.1 SESOI for Three-way Interaction

In the sections above, we derived the SESOI used for the main effect sample-size determination analysis and for the two-way interaction effect-size sensitivity analysis for the binary variables polAff (dem vs. rep) and ewe (TRUE vs. FALSE). subjAttr, however, is a continuous variable that ranges from 1 to 5. Fortunately, in finding a theoretically sound SESOI for the three-way interaction, the same considerations apply as for the main effect and two-way interaction before. We just need to translate these considerations into the continuous metric of subjAttr.

We start by noticing that the complete fixed three-way interaction polAff × ewe × subjAttr is modeled as:

\[ \begin{split} \Delta Duration = \beta_{0} + \\ \beta_{1} \cdot polAff + \beta_{2} \cdot ewe + \beta_{3} \cdot subjAttr + \\ \beta_{4} \cdot (polAff \times ewe) + \beta_{5} \cdot (polAff \times subjAttr) + \beta_{6} \cdot (ewe \times subjAttr) + \\ \beta_{7} \cdot (polAff \times ewe \times subjAttr) \end{split} \]

By rearranging terms, one can show that the two-way interaction polAff × ewe is given by:

\[ Two{-}way\ Interaction_{polAff \times ewe} = \beta_{4} + \beta_{7} \cdot subjAttr \]

Now, let’s calculate this two-way interaction for two individuals who differ in their level of subjective attribution of extreme weather events to climate change. First, an individual who has an average score on subjAttr will show the following two-way interaction effect, with \(\mu_{subjAttr}\) being the sample average of the variable subjAttr:

\[ Effect_{Avg} = \beta_{4} + \beta_{7} \cdot \mu_{subjAttr} \]

Second, we define an individual with a low score on subjAttr as one that shows a subjective attribution of one SD bellow the average. This individual will show the following two-way interaction effect, with \(\sigma_{subjAttr}\) being the SD of subjAttr:

\[ Effect_{Low} = \beta_{4} + \beta_{7} \cdot (\mu_{subjAttr} - \sigma_{subjAttr}) \]

The difference in the two-way interaction effect polAff × ewe between these two individuals is given by:

\[ \begin{split} Effect_{Avg} - Effect_{Low} = \\ \beta_{4} + \beta_{7} \cdot \mu_{subjAttr} - (\beta_{4} + \beta_{7} \cdot (\mu_{subjAttr} - \sigma_{subjAttr})) = \\ \beta_{7} \cdot [\mu_{subjAttr} - (\mu_{subjAttr} - \sigma_{subjAttr})] = \\ \beta_{7} \cdot \sigma_{subjAttr} \end{split} \]

As outlined above in Section Section 3.5.1, we assume that the SESOI for the two-way interaction polAff × ewe is 0.1143 for an average individual (with respect to subjAttr). If the same two-way interaction polAff × ewe shrinks to zero for an individual low in subjAttr, we would consider this interaction effect difference as theoretically relevant (see Figure 14). These assumptions translate to:

\[ Effect_{Avg} - Effect_{Low} = 0.1143 = \beta_{7} \cdot \sigma_{subjAttr} \]

Division by \(\sigma_{subjAttr}\) gives us the SESOI for the three-way interaction in the suitable metric of subjAttr:

\[ SESOI_{polAff \times ewe \times subjAttr} = \frac{0.1143}{\sigma_{subjAttr}} \]

We will assess subjective attribution of extreme weather events to climate change using the same questions, response options, and aggregation as Ogunbode et al. (2019). These authors reported \(\mu_{subjAttr}\) = 3.67 and \(\sigma_{subjAttr}\) = 0.85, resulting in:

\[ SESOI_{polAff \times ewe \times subjAttr} = \frac{0.1143}{0.85} = 0.1345 \]

Figure 14: Assumed ΔDuration values for individuals scoring low (-SD) and average (M) on subjAttr. While individuals scoring average on subjAttr show a two-way interaction effect of polAff × ewe of 0.1143, this effect shrinks to zero for individuals scoring low on subjAttr. These individuals only show the predicted main effect of polAff (whis is 0.1143).

3.6.2 Data Simulation Function

We finally define a function that simulates data for the three-way interaction effect of political affiliation with extreme weather exposure and attribution of extreme weather events to climate change on ΔDuration: FUN_sim_3wayInt. The function will simulate data according to the following model, expressed in lme4-lingo:

deltaDuration ~ polAff * ewe * subjAttr + (1|subj) + (1|trial)

The function FUN_sim_3wayInt takes, among others, the following important arguments (in addition to the arguments discussed for FUN_sim_2wayInt):

• s_mean: Sample mean of subjective attribution. Based on results reported by Ogunbode et al. (2019), we set this value to 3.67.

• s_sd: Sample SD of subjective attribution. Based on results reported by Ogunbode et al. (2019), we set this value to 0.85.

• beta_s: Fixed main effect of subjective attribution. As we are interested in the moderating role of subjective attribution of extreme weather events to climate change, we set this main effect to zero.

• beta_p_s_inx Fixed two-way interaction effect of political affiliation and subjective attribution. In order to accurately model a three-way interaction, one needs to include all two-way interactions in statistical models. Therefore, we include this two-way interaction, but we assume it to be zero.

• beta_e_s_inx Fixed two-way interaction effect of extreme weather exposure and subjective attribution. As before, we include this interaction for accurately modelling the three-way interaction of interest, but we assume this two-way interaction to be zero.

• beta_p_e_s_inx Fixed three-way interaction effect of political affiliation, extreme weather exposure, and subjective attribution. We set this initial value to the SESOI derived above, i.e. 0.1345. We investigate how changing this effect size impacts statistical power, as we are conducting effect-size sensitivity analyses for interaction effects.

The function is defined below:

# define data simulation function
FUN_sim_3wayInt <- function(
  n_subj         =                1000, # number of subjects
  n_subj_prop_p  =           c(.5, .5), # proportion of republican and democrat subjects
  n_subj_prop_e  =           c(.5, .5), # proportion of subjects without and with extreme weather exposure
  n_subj_prop_s  =           c(.5, .5), # proportion of subjects with low and high subjective attribution of EWE to CC
  n_trial        =                  25, # number of trials
  s_mean         =                3.67, # mean of subjAttr, see Ogunbode et al. (2019)
  s_sd           =                0.85, # sd of subjAttr, see Ogunbode et al. (2019)
  beta_0         =                   0, # intercept (grand mean) for deltaDuration
  beta_p         =             0.4/3.5, # main effect of political affiliation (polAff)
  beta_e         =                   0, # main effect of extreme weather exposure (ewe)
  beta_s         =                   0, # main effect of subjective attribution of ewe to climate change (subjAttr)
  beta_p_e_inx   =             0.4/3.5, # two-way interaction effect of polAff and ewe
  beta_p_s_inx   =                   0, # two-way interaction effect of polAff and subjAttr
  beta_e_s_inx   =                   0, # two-way interaction effect of ewe and subjAttr
  beta_p_e_s_inx =      (0.4/3.5)/0.85, # three-way interaction effect of polAff, ewe, and subjAttr
  subj_0         =                 .29, # by-subject random intercept sd for dt carbon
  trial_0        =                 .04, # by-trial random intercept sd
  sigma          =         1*(.29+.04), # residual (error) sd
  
  truncNums      =                TRUE, # should impossible numbers be truncuated?
  setSeed        =                NULL  # seed number to achieve reproducible results. Set to NULL for simulations!
) {
  
  # set seed to achieve reproducible results for demonstration purposes
  set.seed(setSeed)
  
  # simulate data for dwell time on carbon information
  dataSim <- 
    # add random factor subject
    add_random(subj = n_subj) %>% 
    # add random factor trial
    add_random(trial = n_trial) %>% 
    # add between-subject factor political affiliation (with anova contrast)
    add_between("subj", polAff = c("rep", "dem"), .prob = n_subj_prop_p*n_subj, .shuffle = TRUE) %>% 
    add_contrast("polAff", colnames = "X_p", contrast = "anova") %>% 
    # add between-subject factor extreme weather exposure (with anova contrast)
    add_between("subj", ewe = c(FALSE, TRUE), .prob = n_subj_prop_e*n_subj, .shuffle = TRUE) %>% 
    add_contrast("ewe", colnames = "X_e", contrast = "anova") %>% 
    # add between-subject variable subjective attribution of EWE to climate change
    mutate(
      subjAttr = rep(rnorm(n = n_subj, mean = s_mean, sd = s_sd), each = n_trial),
      subjAttr_c = scale(subjAttr, center = TRUE, scale = FALSE)[,1]
    ) %>% 
    # add by-subject random intercept
    add_ranef("subj", S_0 = subj_0) %>% 
    # add by-trial random intercept
    add_ranef("trial", T_0 = trial_0) %>% 
    # add error term
    add_ranef(e_st = sigma) %>% 
    # add response values
    mutate(
      # add together fixed and random effects for each effect
      B_0 = beta_0 + S_0 + T_0,
      B_p = beta_p,
      B_e = beta_e,
      B_s = beta_s,
      B_p_e_inx = beta_p_e_inx,
      B_p_s_inx = beta_p_s_inx,
      B_e_s_inx = beta_e_s_inx,
      B_p_e_s_inx = beta_p_e_s_inx,
      # calculate dv by adding each effect term multiplied by the relevant
      # effect-coded factors and adding the error term
      deltaDuration = 
        B_0 + e_st +
        (X_p * B_p) +
        (X_e * B_e) +
        (subjAttr_c * B_s) +
        (X_p * X_e * B_p_e_inx) +
        (X_p * subjAttr_c * B_p_s_inx) +
        (X_e * subjAttr_c * B_e_s_inx) +
        (X_p * X_e * subjAttr_c * B_p_e_s_inx)
    )
  
  # unset seed
  set.seed(NULL)
  
  # truncuate impossible deltaDurations
  if(truncNums) {
    dataSim <- dataSim %>% 
      mutate(deltaDuration = if_else(deltaDuration < -1, -1,
        if_else(deltaDuration > 1, 1, deltaDuration)))
  }
  
  # run a linear mixed effects model and check summary
  mod <- lmer(
    deltaDuration ~ polAff*ewe*subjAttr_c + (1 | subj) + (1 | trial),
    data = dataSim
  )
  mod.sum <- summary(mod)

  # get results in tidy format
  mod.broom <- broom.mixed::tidy(mod)

  return(list(
    dataSim = dataSim,
    modelLmer = mod,
    modelResults = mod.broom
  ))
  
}

We call the function once and extract the results of this single simulation:

# Note to myself: Consider setting beta_p = 0.4/3.5, beta_p_e_inx = 0.4/3.5,
# and beta_p_e_s_inx to 2 * 0.4/3.5/(2*.85).
# This way, the interaction effect polAff:ewe is 0.4/3.5 for individuals
# with an average subjAttr (subjAttr_c = 0). For individuals with
# subjAttr = mean - SD, polAff:ewe is 0. For individuals with subjAttr = mean + SD,
# polAff:ewe is  2 * 0.4/3.5.

out <- FUN_sim_3wayInt(
  n_subj         =                1000, # number of subjects
  n_subj_prop_p  =           c(.5, .5), # proportion of republican and democrat subjects
  n_subj_prop_e  =           c(.5, .5), # proportion of subjects without and with extreme weather exposure
  n_subj_prop_s  =           c(.5, .5), # proportion of subjects with low and high subjective attribution of EWE to CC
  n_trial        =                  25, # number of trials
  s_mean         =                3.67, # mean of subjAttr, see Ogunbode et al. (2019)
  s_sd           =                0.85, # sd of subjAttr, see Ogunbode et al. (2019)
  beta_0         =                   0, # intercept (grand mean) for deltaDuration
  beta_p         =             0.4/3.5, # main effect of political affiliation (polAff)
  beta_e         =                   0, # main effect of extreme weather exposure (ewe)
  beta_s         =                   0, # main effect of subjective attribution of ewe to climate change (subjAttr)
  beta_p_e_inx   =             0.4/3.5, # two-way interaction effect of polAff and ewe
  beta_p_s_inx   =                   0, # two-way interaction effect of polAff and subjAttr
  beta_e_s_inx   =                   0, # two-way interaction effect of ewe and subjAttr
  beta_p_e_s_inx =      (0.4/3.5)/0.85, # three-way interaction effect of polAff, ewe, and subjAttr
  subj_0         =                 .29, # by-subject random intercept sd for dt carbon
  trial_0        =                 .04, # by-trial random intercept sd
  sigma          =         1*(.29+.04), # residual (error) sd
  
  truncNums      =                TRUE, # should impossible numbers be truncuated?
  setSeed        =                 123  # seed number to achieve reproducible results. Set to NULL for simulations!
)

# Get results table
resultsTable <- out$modelResults %>% 
  select(-c(std.error, statistic, df)) %>% 
  mutate(across(where(is_double), ~ round(.x, 4))) %>% 
  knitr::kable()
formulaUsedForFit <- paste(as.character(formula(out$modelLmer))[c(2,1,3)], collapse = " ")

# Create predictions plot

# refit model to dispaly subjAttr levels in original metric (not mean centered)
m <- lmer(
  deltaDuration ~ polAff*ewe*subjAttr + (1 | subj) + (1 | trial),
  data = out$dataSim
)
# define spotlights for spotlight analysis
spotlights <- c(3.67 - 0.85, 3.67, 3.67 + 0.85)

# create plot showing predictions
p.demo.3wayInt.pred <- predict_response(m, terms = c("ewe", "polAff", "subjAttr[spotlights]"))
p.demo.3wayInt <- p.demo.3wayInt.pred %>% 
  plot(colors = c("red", "dodgerblue")) +
  coord_cartesian(ylim = c(-.25, .25)) +
  theme_bw()

jpeg(
  file = "../images/pa_demo3wayInt.jpeg",
  width = 9.1, height = 6.3, units = "in", res = 600
)
print(p.demo.3wayInt)
invisible(dev.off())

Figure 15 visualizes predictions based on this single simulation and Table 5 summarizes the statistical results of fitting the actual model used in data generation to the simulated data.

Figure 15: Visual representation of results of one simulation created using FUN_sim_3wayInt. Points indicate the predicted means surrounded by 95% Confidence Intervals. Panels indicate predictions for different values of subjAttr (Mean - SD, Mean, and Mean + SD).

Table 5: Statistical results of one simulation created using FUN_sim_3wayInt. Data was fit using deltaDuration ~ polAff * ewe * subjAttr_c + (1 | subj) + (1 | trial). Note that subjAttr is mean centered to ease interpretation of lower-level interactions and main effects.

effect	group	term	estimate	p.value
fixed	NA	(Intercept)	-0.0022	0.8478
fixed	NA	polAff.dem-rep	0.1149	0.0000
fixed	NA	ewe.TRUE-FALSE	-0.0023	0.9031
fixed	NA	subjAttr_c	0.0091	0.4141
fixed	NA	polAff.dem-rep:ewe.TRUE-FALSE	0.1426	0.0001
fixed	NA	polAff.dem-rep:subjAttr_c	0.0187	0.4025
fixed	NA	ewe.TRUE-FALSE:subjAttr_c	0.0142	0.5239
fixed	NA	polAff.dem-rep:ewe.TRUE-FALSE:subjAttr_c	0.1111	0.0130
ran_pars	subj	sd__(Intercept)	0.2849	NA
ran_pars	trial	sd__(Intercept)	0.0323	NA
ran_pars	Residual	sd__Observation	0.3240	NA

3.6.3 Power Simulation

In the following code, the simulations are calculated. We do not recommend executing this code junk as it takes several hours to run.

FUN_sim_3wayInt_pwr <- function(sim, ...){
  out <- FUN_sim_3wayInt(...)
  modelResults <- out$modelResults %>% 
    mutate(sim = sim) %>% 
    relocate(sim)
  return(modelResults)
}

# How many simulations should be run?
n_sims <- 1000

# What are the breaks for number of subjects we would like to calculate power for?
breaks_subj <- c(900, 950, 1000)

# What are the breaks for SESOI?
breaks_sesoi <- (0.4/3.5)/0.85 * seq(1, 2, .25)

# What are the breaks for different error SDs?
breaks_sigma <- c((.29+.04), 2*(.29+.04))

res_3wayInt <- tibble()
for (s in seq_along(breaks_sigma)) {
  
  res_sesoi <- tibble()
  for (sesoi in seq_along(breaks_sesoi)) {
    
    res_nSubj <- tibble()
    for (nSubj in seq_along(breaks_subj)) {
      
      # Give feedback regarding which model is simulated
      cat(paste0(
        "Simulation:\n",
        "  sigma = ", round(breaks_sigma[s], 4), "\n",
        "  sesoi = ", round(breaks_sesoi[sesoi], 4), "\n",
        "  nSubject = ", breaks_subj[nSubj], "\n"
      ))
      
      # Start timer
      cat(paste0("Start date time: ", lubridate::now(), "\n"))
      tic()
      
      # Loop over simulations
      pwr <- map_df(
        1:n_sims, 
        FUN_sim_3wayInt_pwr,
        n_subj = breaks_subj[nSubj],
        beta_p_e_s_inx = breaks_sesoi[sesoi],
        sigma = breaks_sigma[s]
      )
      
      # Stop timer and calculate elapsed time
      elapsed_time <- toc(quiet = TRUE)
      elapsed_seconds <- elapsed_time$toc - elapsed_time$tic
      elapsed_minutes <- elapsed_seconds / 60
      cat(paste0("End date time: ", lubridate::now(), "\n"))
      cat("Elapsed time: ", elapsed_minutes, " minutes\n\n")
      
      # Add number of subjects to pwr
      pwr <- pwr %>% 
        mutate(
          nSubjects = breaks_subj[nSubj],
          sesoi = breaks_sesoi[sesoi],
          sigma = breaks_sigma[s]
        )
      
      # Add results to the results table
      res_nSubj <- res_nSubj %>%
        rbind(pwr)
    }
    
    # Add results to the results table
    res_sesoi <- res_sesoi %>% 
      rbind(res_nSubj)
  }
  
  # Add results to the results table
  res_3wayInt <- res_3wayInt %>% 
    rbind(res_sesoi)
  
}

res_3wayInt.summary <- res_3wayInt %>% 
  filter(term == "polAff.dem-rep:ewe.TRUE-FALSE:subjAttr_c") %>% 
  group_by(sigma, sesoi, nSubjects) %>% 
  summarise(
    power = mean(p.value < 0.05),
    ci.lower = binom.confint(power*n_sims, n_sims, methods = "exact")$lower,
    ci.upper = binom.confint(power*n_sims, n_sims, methods = "exact")$upper,
    .groups = 'drop'
  ) %>% 
  mutate(
    sigma_fact = factor(format(round(sigma, 4), nsmall = 4)),
    sigma_level = match(sigma_fact, levels(sigma_fact)),
    sesoi_fact = factor(format(round(sesoi, 4), nsmall = 4)),
    sesoi_level = match(sesoi_fact, levels(sesoi_fact))
  )

# Save results in a list object
time <- format(Sys.time(), "%Y%m%d_%H%M")
fileName <- paste0("res_3wayInt", "_", time, ".RDS")
saveRDS(
  list(
    res_3wayInt = res_3wayInt,
    res_3wayInt.summary = res_3wayInt.summary
  ),
  file = file.path("../powerSimulationsOutput", fileName)
)

We retrieve pre-run results:

# Load power simulation data
resList_3wayInt <- readRDS(file.path("../powerSimulationsOutput", "res_3wayInt_20240814_1612.RDS")) 
resList_3wayInt.summary <- resList_3wayInt$res_3wayInt.summary

# Extract power values for some specific effect sizes at N = 1000
powerValues <- resList_3wayInt.summary %>% 
  filter(sigma_fact == "0.3300") %>% 
  filter(sesoi_fact == "0.1345") %>% 
  filter(nSubjects == 950) %>% 
  mutate(power_str = paste0(round(power*100, 2), "%")) %>% 
  pull(power_str)

# Extract number of simulations
label_nSimulations <- resList_3wayInt$res_3wayInt$sim %>% n_distinct()

# Repeat breaks_sesoi
breaks_sesoi <- (0.4/3.5)/(0.85) * seq(1, 2, .25)

Figure 16 displays the distribution of estimated fixed effects across all simulations. The figure shows that the estimated fixed effects are close to the true ones provided as input in the data simulation function, validating that simulations worked as expected.

Figure 16: Distribution of estimated fixed effects resulting from 1000 simulations for the model deltaDuration ~ polAff * ewe * subjAttr_c + (1 | subj) + (1 | trial). Shaded area represent densities, annotated points indicate medians, and thick and thin lines represent 66% and 95% quantiles.

Figure 17 shows results of our effect-size sensitivity analyses. We plot statistical power (y-axis) for different effect sizes (x-axis), taking into account different assumptions for the error SD (color) and sample size (panel). Regarding the latter, we report results not only for the full sample size we aim for (N = 1000), but also for sample sizes taking into account different participant exclusion-rates due to exclusion criteria defined in the Registered Report.

Figure 17: Power curves for the three-way interaction polAff × ewe × subjAttr. Points represent simulated power surrounded by a 95%-CI based on 1000 simulations with α = 0.05. Note that, in contrast to Figure 9, the x-axsis represents different effect sizes, starting from the defined SESOI, while the panels represent different sample sizes, taking into account participant exclusion-rates of 10% (N = 900), 5% (N = 950), and 0% (N = 1000). Note that estimates are displayed with a slight shift along the x-axis to reduce overlap. Since we only consider effect sizes equal to or greater than the SESOI to be of relevance, power curves were only examined for positive effect sizes, which is why only positive values are shown.

# Get power values of interest
chosenN <- 950
chosenSigma <- c("0.3300", "0.6600")
chosenSESOI <- c("0.1681", "0.2017")
powerValues <- resList_3wayInt.summary %>% 
  filter(sigma_fact %in% chosenSigma) %>% 
  filter(sesoi_fact %in% chosenSESOI) %>% 
  filter(nSubjects %in% chosenN)

# Get lower CI for power values for more liberal and more conservative sigma assumptions
powerValues_sigmaLib <- powerValues %>% 
  filter(sigma_fact == "0.3300")
powerValues_sigmaCons <- powerValues %>% 
  filter(sigma_fact == "0.6600")

# Interpolate the effect sizes at which we achieve 95% power
eff_sigmaLib <- approx(powerValues_sigmaLib$ci.lower, powerValues_sigmaLib$sesoi, xout = 0.95)$y
eff_sigmaCons <- approx(powerValues_sigmaCons$ci.lower, powerValues_sigmaCons$sesoi, xout = 0.95)$y
# Round results for display in text
eff_sigmaLib_txt_3way <- format(round(eff_sigmaLib, 4), nsmall = 4)
eff_sigmaCons_txt_3way <- format(round(eff_sigmaCons, 4), nsmall = 4)

To assess the smallest effect size that can be detected with 95% statistical power, we inspect the lower bounds of the 95%-CI power estimates in Figure 17. Specifically, we focus on the power simulation results for N = 950, which takes into account a participant exclusion rate of 5%. There, we interpolate between the two point estimates that lie just below and above the 95% power line, i.e., between the power estimates for effect sizes 0.1681 and 0.2017. Assuming an error SD of 0.3300, we achieve 95% statistical power to detect a three-way interaction effect of at least 0.1728. For a more conservative error SD of 0.6600, this smallest detectable effect size is only marginally higher (0.1890).

3.6.3.1 Increased Sample Size

We additionally conduct an effect-size sensitivity analysis for increased sample sizes.

FUN_sim_3wayInt_pwr <- function(sim, ...){
  out <- FUN_sim_3wayInt(...)
  modelResults <- out$modelResults %>% 
    mutate(sim = sim) %>% 
    relocate(sim)
  return(modelResults)
}

# How many simulations should be run?
n_sims <- 1000

# What are the breaks for number of subjects we would like to calculate power for?
breaks_subj <- c(1500, 1750, 2000)

# What are the breaks for SESOI?
breaks_sesoi <- (0.4/3.5)/0.85 * seq(1, 2, .25)

# What are the breaks for different error SDs?
breaks_sigma <- c((.29+.04), 2*(.29+.04))

res_3wayInt <- tibble()
for (s in seq_along(breaks_sigma)) {
  
  res_sesoi <- tibble()
  for (sesoi in seq_along(breaks_sesoi)) {
    
    res_nSubj <- tibble()
    for (nSubj in seq_along(breaks_subj)) {
      
      # Give feedback regarding which model is simulated
      cat(paste0(
        "Simulation:\n",
        "  sigma = ", round(breaks_sigma[s], 4), "\n",
        "  sesoi = ", round(breaks_sesoi[sesoi], 4), "\n",
        "  nSubject = ", breaks_subj[nSubj], "\n"
      ))
      
      # Start timer
      cat(paste0("Start date time: ", lubridate::now(), "\n"))
      tic()
      
      # Loop over simulations
      pwr <- map_df(
        1:n_sims, 
        FUN_sim_3wayInt_pwr,
        n_subj = breaks_subj[nSubj],
        beta_p_e_s_inx = breaks_sesoi[sesoi],
        sigma = breaks_sigma[s]
      )
      
      # Stop timer and calculate elapsed time
      elapsed_time <- toc(quiet = TRUE)
      elapsed_seconds <- elapsed_time$toc - elapsed_time$tic
      elapsed_minutes <- elapsed_seconds / 60
      cat(paste0("End date time: ", lubridate::now(), "\n"))
      cat("Elapsed time: ", elapsed_minutes, " minutes\n\n")
      
      # Add number of subjects to pwr
      pwr <- pwr %>% 
        mutate(
          nSubjects = breaks_subj[nSubj],
          sesoi = breaks_sesoi[sesoi],
          sigma = breaks_sigma[s]
        )
      
      # Add results to the results table
      res_nSubj <- res_nSubj %>%
        rbind(pwr)
    }
    
    # Add results to the results table
    res_sesoi <- res_sesoi %>% 
      rbind(res_nSubj)
  }
  
  # Add results to the results table
  res_3wayInt <- res_3wayInt %>% 
    rbind(res_sesoi)
  
}

res_3wayInt.summary <- res_3wayInt %>% 
  filter(term == "polAffX_p:eweX_e:subjAttr_c") %>% 
  group_by(sigma, sesoi, nSubjects) %>% 
  summarise(
    power = mean(p.value < 0.05),
    ci.lower = binom.confint(power*n_sims, n_sims, methods = "exact")$lower,
    ci.upper = binom.confint(power*n_sims, n_sims, methods = "exact")$upper,
    .groups = 'drop'
  ) %>% 
  mutate(
    sigma_fact = factor(format(round(sigma, 4), nsmall = 4)),
    sigma_level = match(sigma_fact, levels(sigma_fact)),
    sesoi_fact = factor(format(round(sesoi, 4), nsmall = 4)),
    sesoi_level = match(sesoi_fact, levels(sesoi_fact))
  )

# Save results in a list object
time <- format(Sys.time(), "%Y%m%d_%H%M")
fileName <- paste0("res_3wayInt_increasedN", "_", time, ".RDS")
saveRDS(
  list(
    res_3wayInt = res_3wayInt,
    res_3wayInt.summary = res_3wayInt.summary
  ),
  file = file.path("../powerSimulationsOutput", fileName)
)

# Load power simulation data
resList_3wayInt <- readRDS(file.path("../powerSimulationsOutput", "res_3wayInt_increasedN_20250308_1649.RDS"))
resList_3wayInt.summary <- resList_3wayInt$res_3wayInt.summary

# Extract power values for some specific effect sizes at N = 1000
powerValues <- resList_3wayInt.summary %>% 
  filter(sigma_fact == "0.3300") %>% 
  filter(sesoi_fact == "0.1345") %>% 
  filter(nSubjects == 950) %>% 
  mutate(power_str = paste0(round(power*100, 2), "%")) %>% 
  pull(power_str)

# Extract number of simulations
label_nSimulations <- resList_3wayInt$res_3wayInt$sim %>% n_distinct()

# Repeat breaks_sesoi
breaks_sesoi <- (0.4/3.5)/(0.85) * seq(1, 2, .25)

Figure 18: Power curves for the three-way interaction polAff × ewe × subjAttr. Points represent simulated power surrounded by a 95%-CI based on 1000 simulations with α = 0.05. Note that, in contrast to Figure 9, the x-axsis represents different effect sizes, starting from the defined SESOI, while the panels represent different sample sizes. Note that estimates are displayed with a slight shift along the x-axis to reduce overlap. Since we only consider effect sizes equal to or greater than the SESOI to be of relevance, power curves were only examined for positive effect sizes, which is why only positive values are shown.

3.7 Conclusion

Using power simulations with mixed-effects models for sample-size determination and effect-size sensitivity analyses, we show that with a final sample size of N = 950:

1. We will achieve \(\ge\) 95% statistical power to detect a main effect of political affiliation.

2. We will be able to detect a two-way interaction effect of political affiliation with extreme weather exposure of at least 0.1619 with \(\ge\) 95% statistical power.

3. We will be able to detect a three-way interaction effect of political affiliation with extreme weather exposure and subjective attribution of extreme weather events to climate change of at least 0.1890 with \(\ge\) 95% statistical power.

By increasing the planned sample size to N = 2’000, we will be able to detect all of the defined smallest effect sizes of interest of the hypothesized main effect and the two- and three-way interaction effects with a statistical power of at least 95%.

Expand for Session Info

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.4.0 (2024-04-24)
 os       macOS Sonoma 14.4.1
 system   aarch64, darwin20
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       Europe/Zurich
 date     2024-09-19
 pandoc   3.1.11 @ /usr/local/bin/ (via rmarkdown)
 quarto   1.5.45 @ /usr/local/bin/quarto

─ Packages ───────────────────────────────────────────────────────────────────
 package     * version  date (UTC) lib source
 binom       * 1.1-1.1  2022-05-02 [1] CRAN (R 4.4.0)
 dplyr       * 1.1.4    2023-11-17 [1] CRAN (R 4.4.0)
 faux        * 1.2.1    2023-04-20 [1] CRAN (R 4.4.0)
 forcats     * 1.0.0    2023-01-29 [1] CRAN (R 4.4.0)
 fs          * 1.6.4    2024-04-25 [1] CRAN (R 4.4.0)
 ggdist      * 3.3.2    2024-03-05 [1] CRAN (R 4.4.0)
 ggeffects   * 1.7.0    2024-06-20 [1] CRAN (R 4.4.0)
 ggplot2     * 3.5.1    2024-04-23 [1] CRAN (R 4.4.0)
 ggpubr      * 0.6.0    2023-02-10 [1] CRAN (R 4.4.0)
 ggrepel     * 0.9.5    2024-01-10 [1] CRAN (R 4.4.0)
 ggthemes    * 5.1.0    2024-02-10 [1] CRAN (R 4.4.0)
 lme4        * 1.1-35.5 2024-07-03 [1] CRAN (R 4.4.0)
 lmerTest    * 3.1-3    2020-10-23 [1] CRAN (R 4.4.0)
 lubridate   * 1.9.3    2023-09-27 [1] CRAN (R 4.4.0)
 Matrix      * 1.7-0    2024-03-22 [1] CRAN (R 4.4.0)
 purrr       * 1.0.2    2023-08-10 [1] CRAN (R 4.4.0)
 readr       * 2.1.5    2024-01-10 [1] CRAN (R 4.4.0)
 sessioninfo * 1.2.2    2021-12-06 [1] CRAN (R 4.4.0)
 stringr     * 1.5.1    2023-11-14 [1] CRAN (R 4.4.0)
 tibble      * 3.2.1    2023-03-20 [1] CRAN (R 4.4.0)
 tictoc      * 1.2.1    2024-03-18 [1] CRAN (R 4.4.0)
 tidyr       * 1.3.1    2024-01-24 [1] CRAN (R 4.4.0)
 tidyverse   * 2.0.0    2023-02-22 [1] CRAN (R 4.4.0)
 usmap       * 0.7.1    2024-03-21 [1] CRAN (R 4.4.0)

 [1] /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library

──────────────────────────────────────────────────────────────────────────────

References

Berger, Sebastian, and Annika M Wyss. 2021. “Measuring Pro-Environmental Behavior Using the Carbon Emission Task.” Journal of Environmental Psychology 75: 101613. https://doi.org/https://doi.org/10.1016/j.jenvp.2021.101613.

Giner-Sorolla, Roger, Amanda K. Montoya, Alan Reifman, Tom Carpenter, Neil A. Lewis, Christopher L. Aberson, Dries H. Bostyn, et al. 2024. “Power to Detect What? Considerations for Planning and Evaluating Sample Size.” Personality and Social Psychology Review, February, 10888683241228328. https://doi.org/10.1177/10888683241228328.

IPCC, ed. 2023. “Weather and Climate Extreme Events in a Changing Climate.” In, 1513–1766. Cambridge: Cambridge University Press. https://doi.org/10.1017/9781009157896.013.

Konisky, David M., Llewelyn Hughes, and Charles H. Kaylor. 2016. “Extreme Weather Events and Climate Change Concern.” Climatic Change 134 (4): 533–47. https://doi.org/10.1007/s10584-015-1555-3.

Ogunbode, Charles A., Christina Demski, Stuart B. Capstick, and Robert G. Sposato. 2019. “Attribution Matters: Revisiting the Link Between Extreme Weather Experience and Climate Change Mitigation Responses.” Global Environmental Change 54 (January): 31–39. https://doi.org/10.1016/j.gloenvcha.2018.11.005.

Reeck, Crystal, Daniel Wall, and Eric J. Johnson. 2017. “Search Predicts and Changes Patience in Intertemporal Choice.” Proceedings of the National Academy of Sciences 114 (45): 11890–95. https://doi.org/10.1073/pnas.1707040114.

Salthouse, Timothy A. 2000. “Aging and Measures of Processing Speed.” Biological Psychology 54 (1): 35–54. https://doi.org/10.1016/S0301-0511(00)00052-1.

Willemsen, Martijn C., and Eric J. Johnson. 2019. “(Re)visiting the Decision Factory: Obswerving Cognition with MouselabWEB.” In, edited by Michael Schulte-Mecklenbeck, Anton Kuehberger, and Joseph G. Johnson, 2nd ed., 76–95. New York: Routledge. https://doi.org/10.4324/9781315160559.